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Chapter  1 
Abstract 


During  the  period  from  April  2003  to  March  2006,  this  research  had  been  progressing 
well  as  planned  toward  the  ultimate  goal  of  simulating  the  mistuned  rotor  with  fully- 
coupled  fluid  structure  interaction.  This  is  a  multidisciplinary  comprehensive  project 
that  needs  components  from  both  fluid  and  structure  dynamics.  The  following  im¬ 
portant  capabilities  have  been  achieved  and  the  detailed  results  are  presented  in  the 
report. 

1)  A  new  accurate  and  efficient  Riemann  solver,  Zha  E-CUSP2  scheme,  has  been 
developed  and  applied  to  moving  grid  systems  for  fully  coupled  fluid-structure  inter¬ 
action.  The  schemes  are  validated  to  possess  very  low  numerical  diffusion  and  be  able 
to  capture  crisp  shock  and  contact  discontinuities.  The  robustness  and  efficiency  of 
the  scheme  is  essential  for  the  calculation  of  fluid-structural  interactions. 

2)  The  fully  coupled  fluid-structural  interaction  methodology  has  been  developed 
and  is  successfully  applied  to  predict  the  2D  and  3D  transonic  wing  flutter.  For 
the  fluid-structural  interaction,  an  implicit  time  marching  method  with  dual  time 
stepping  algorithm  and  unfactored  Gauss-Seidel  line  relaxation  is  employed  to  achieve 
fast  convergence  rate.  For  the  2D  cases,  the  exact  structural  equations  are  solved. 
For  the  3D  AGARD  wing,  a  modal  approach  is  used  to  be  consistent  with  the  Subsets 
of  Nominal  Modes  (SNM)  model  of  the  Mistuned  rotor. 

3)  A  non-reflective  boundary  condition  based  on  Navier-Stokes  equations  in  gen¬ 
eralized  coordinates  has  been  developed  to  accurately  treat  the  boundaries  of  the 
unsteady  flows  to  remove  the  reflective  waves. 


4)  The  transient  response  (time  domain)  structural  vibration  model  for  mistuned 
rotor  bladed  disk  based  on  the  efficient  SNM  model  has  been  developed.  The  vi¬ 
bration  response  results  predicted  by  the  SNM  model  for  a  full  annulus  bladed  disk 
with  blade  frequency  variation  agree  very  well  with  the  results  predicted  by  the  finite 
element  model. 


3 


4 


CHAPTER  1.  ABSTRACT 


5)  The  code  has  been  intensively  validated  with  several  eases  including  steady  state 
2D  transonic  airfoil  and  3D  wing,  unsteady  vortex  shedding  of  a  stationary  cylinder, 
induced  vibration  of  a  cylinder,  forced  vibration  of  a  pitching  airfoil,  induced  vibration 
and  flutter  boundary  of  2D  NACA  64A010  transonic  airfoil,  3D  plate  wing  structural 
response.  The  predicted  results  agree  well  with  benchmark  experimental  results  or 
the  results  calculated  by  a  finite  element  solver  for  structural  response.  The  limited 
cycle  oscillation  (LCO)  is  captured. 

6)  The  full  3D  AGARD  wing  flutter  boundary  is  calculated  and  agree  well  with 
the  experiment.  The  “sonic  dip”  phenomenon  is  captured. 

This  solver  based  on  the  fully  coupled  fluid-structural  interaction  is  ready  to  calculate 
the  mistuned  rotor  flutter  and  forced  response.  However,  since  the  funding  is  only  3 
years,  which  is  one  year  shorter  than  the  4  years  time  period  originally  proposed,  the 
mistuned  rotor  simulation  is  not  finished  and  will  be  completed  in  future  when  the 
funding  is  available. 

In  this  research  project,  we  have  5  journal  papers,  2  papers  submitted  for  journal 
publications,  and  12  conference  papers. 
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Chapter  2 
Introduction 


Flow  induced  structural  vibration  is  one  of  the  most  critical  technical  problems  affect¬ 
ing  the  readiness  of  the  US  Air  Force  fleet  today.  Due  to  the  extremely  complicated 
non-linear  flow-structure  interaction  phenomena,  there  is  a  lack  of  high  fidelity  com¬ 
putational  tools  to  study  the  basic  physics  and  to  predict  the  structural  failure.  The 
problems  are  more  complex  for  the  aircraft  engine  turbomachinery  than  airframe  be¬ 
cause  the  turbomachinery  has  numerous  blades  and  is  much  more  complicated  than 
the  one  or  two  wings  of  the  external  airframe. 

For  aircraft  engine  turbomachinery,  the  most  general  flow  induced  vibration  prob¬ 
lem  is  the  mistuning  problem.  The  mistuning  refers  to  the  fact  that,  in  a  bladed  disk, 
blade  responses  can  significantly  differ  from  each  other  due  to  small  geometric  varia¬ 
tions  from  blade  to  blade.  The  small  variations  typically  result  from  within-tolerance 
manufacturing  imperfections  or  in-service  wear-offs,  which  are  difficult  to  eliminate  in 
the  production  process  or  during  the  life  span  of  an  engine.  In  other  words,  virtually 
all  turbomachinery  bladed  disks  are  mistuned.  However,  up  to  date,  almost  all  of  the 
fluid-structure  interaction  research  work  focuses  on  the  tunned  systems  because  the 
mistuned  systems  are  much  more  challenging. 

This  research  is  to  develop  a  methodology  to  couple  a  state  of  the  art  CFD  code, 
RANS3D  (3D  unsteady  Reynolds  averaged  NS  solver),  with  a  state  of  the  art  struc¬ 
tural  model,  SNM  (Subset  of  Nominal  Modes),  for  a  mistuned  bladed  disk.  The 
aerodynamic  force  and  blade  motion  of  a  full  annulus  rotor  are  unknown  variables 
and  are  simultaneously  solved  within  each  time  step.  No  prescribed  blade  motion  is 
used  to  accurately  represent  the  coupled  system.  The  process  proceeds  in  temporal 
direction  step  by  step  until  it  reaches  desired  solutions  for  forced  response  or  flutter. 

The  methodology  developed  in  this  research  is  based  on  the  strategy  of  computing 
frameworks  recently  suggested  by  Melville[7],  The  computing  frameworks  is  to  fully 
take  advantages  of  the  state  of  the  art  of  each  individual  discipline  and  the  multi¬ 
disciplinary  computation  is  connected  through  standard  interfaces.  In  this  research, 
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the  multi-disciplines  involved  are  the  CFD  model  and  the  structural  model.  Since 
the  methodologies  of  each  individual  discipline  have  been  developed  and  matured 
independently,  each  methodology  will  usually  use  different  grid  generation  and  time 
marching  scheme  to  achieve  the  highest  possible  accuracy  and  efficiency.  A  com¬ 
mon  interface  to  conserve  the  energy  between  the  fluid  aerodynamic  forcing  and  the 
structural  deformation  is  necessary  to  combine  the  CFD  and  structural  solvers.  The 
strategy  of  the  computing  frameworks  will  allow  the  two  disciplines  to  continue  to 
develop  their  state  of  the  art  and  ensure  that  the  multi-disciplinary  solver  will  always 
be  at  the  forefront  of  the  technology. 

Currently,  the  method  used  to  predict  the  aerodynamic  damping  of  mistuned 
bladed  disks  is  based  on  the  Influence  Coefficient  Technique  which  requires  prescribed 
blade  motion  (Silkowski  et.  al.,  2001)  [8].  This  technique  assumes  no  structural  cou¬ 
pling  between  blades  and  does  not  account  for  the  feedback  from  the  mistuned  struc¬ 
ture  to  the  fluid  field  or  vice  versa.  The  accuracy  of  such  approach  is  questionable 
when  the  structural  coupling  is  important  or  highly  nonlinear  flow  conditions  exist. 
For  instance,  at  stall  flutter  region,  there  may  exist  large  oscillating  separation,  ro¬ 
tating  stall  cells,  shock  motion,  oscillating  tip  vortex,  etc.  Under  this  circumstance, 
it  is  not  feasible  to  prescribe  the  blade  motion  caused  by  the  complicated  interaction 
between  the  flow  and  the  structure. 

To  simulate  the  vibratory  response  of  a  mistuned  bladed  disk,  the  sector  model  used 
for  a  tuned  system  is  no  longer  applicable.  A  model  that  simulates  the  whole  bladed 
disk  is  needed  to  take  into  account  the  asymmetric  pattern  of  the  blade  vibration 
around  the  full  annulus  (Hilbert  and  Blair,  2001)  [9],  which  includes  the  differences 
in  amplitudes,  unequal  phase  differences  between  adjacent  blades,  and  sometimes 
changes  in  blade  mode  shapes. 

There  are  increasing  efforts  recently  to  couple  a  computational  fluid  dynamics(CFD) 
solver  with  a  structural  solver,  in  particular  for  external  airframe  problems.  Bendik- 
sen  et  al.  pioneered  the  research  by  using  an  explicit  CFD  code  coupled  with  a 
structural  integrator  based  on  the  convolution  integral  to  obtain  the  flutter  boundary 
for  a  NACA  64A010  air  foil  [10],  Alonso  and  Jameson  etc.  developed  a  dual-time  step 
coupled  aeroelastic  solver  for  2D  airfoilsfll,  12],  Similar  techniques  was  applied  to 
3D  wing  with  an  finite  element  structural  model  by  Liu  etc[13].  Melville  et  al.  devel¬ 
oped  a  fully  implicit  method  based  on  Beam- Warming  scheme  with  a  coupled  modal 
structural  solver  for  a  3D  wing[14]. 

Due  to  the  difficulties  in  turbomachinery,  calculation  of  fluid  induced  vibration 
based  on  fully  coupled  fluid-structure  interaction  has  began  only  very  recently.  The 
research  group  in  UK  led  by  Dr.  Imregun  has  made  notable  progresses  (Breard,  et 
al.  1999,  Sayma  2001a,  Sayma  2001b,  Sayma  2001c,  Sayma  2001d)  [15] [16] [17]  [18]. 
They  carried  out  full  annulus  and  multiblade  row  computations  for  forced  response 
and  flutter.  The  flow  solutions  are  primarily  based  on  Euler  equations  with  limited 
3D  Navier-Stokes  modeling  and  the  structural  response  is  simulated  by  mode  shapes 
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of  bladed  disks.  Nonlinear  friction  constraints  are  considered  in  their  study;  however, 
no  mistuning  study  was  attempted.  In  addition,  in  their  work,  a  mode  superposi¬ 
tion  of  the  structure  is  incorporated  into  a  finite  element  CFD  solver.  This  method 
hence  may  not  best  fit  the  computing  frameworks [7].  In  2002,  Doi  and  Alonso[19] 
applied  their  dual  time  stepping  CFD  algorithm[ll,  12]  with  the  structural  solver  of 
MSC/NASTRAN  to  simulate  a  tuned  compressor  rotor  fluid-structural  interaction. 
In  the  model  of  Doi  and  Alonso  [19],  the  fluid  and  structural  models  are  closely  cou¬ 
pled  but  structural  deformation  is  lagged.  This  method  may  be  limited  to  first-order 
accuracy  in  time  regardless  of  the  temporal  accuracy  of  the  individual  solvers  [14]. 

To  achieve  the  research  goal,  the  numerical  strategy  is  given  in  the  next  chapter. 
The  details  of  the  numerical  algorithms  are  given  in  chapter  4. 


CHAPTER  2.  INTRODUCTION 


Chapter  3 

Numerical  Strategy 


The  challenges  for  high  cycle  fatigue  prediction  based  on  a  fully  coupled  fluid-structural 
interaction  are  two-folds,  efficiency  and  accuracy.  The  fully  coupled  fluid-structural 
interaction  needs  to  iterate  the  flow  and  structural  solver  within  each  physical  time 
step.  It  is  a  very  CPU  time  consuming  task  and  the  dynamic  response  of  the  system 
is  sensitive  to  the  numerical  dissipation  introduced  by  the  numerical  scheme.  Conse¬ 
quently,  it’s  required  that  the  numerical  scheme  is  able  to  model  the  flow  field  with 
high  efficiency  and  low  numerical  diffusion. 


3.1  Low  Diffusion  High  Efficiency  Upwind  Scheme 

Recently,  there  have  been  many  efforts  to  develop  efficient  Riemann  solvers  using 
scalar  dissipation  instead  of  matrix  dissipation.  For  the  scalar  dissipation  Riemann 
solver  schemes,  there  are  generally  two  types:  H-CUSP  schemes  and  E-CUSP  schemes 
[20,  21,  22].  The  abbreviation  CUSP  stands  for  “convective  upwind  and  split  pres¬ 
sure”,  named  by  [20,  21,  22].  The  H-CUSP  schemes  have  the  total  enthalpy  from  the 
energy  equation  in  their  convective  vector,  while  the  E-CUSP  schemes  use  the  total 
energy  in  the  convective  vector.  Liou’s  AUSM  family  schemes  [23],  Van  Leer-Hanel 
scheme  [24],  and  Edwards’s  LDFSS  schemes  [25,  26]  belong  to  the  H-CUSP  group. 

The  H-CUSP  schemes  may  have  the  advantage  of  better  conserving  the  total  en¬ 
thalpy  for  steady  state  flows.  However,  from  the  characteristic  theory  point  of  view, 
the  H-CUSP  schemes  are  not  fully  consistent  with  the  disturbance  propagation  direc¬ 
tions,  which  may  affect  the  stability  and  robustness  of  the  schemes  [1].  The  H-CUSP 
scheme  may  have  more  inconsistencies  when  it  is  extended  to  the  moving  grid  system. 
It  will  leave  the  pressure  term  multiplied  by  the  grid  velocity  in  the  energy  flux,  which 
cannot  be  contained  in  the  total  enthalpy,  and  must  therefore  be  treated  as  part  of 
the  pressure  term.  From  a  characteristic  point  of  view,  it  is  not  obvious  how  to  treat 
this  term  in  a  consistent  manner. 
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In  this  research,  we  developed  an  efficient  E-CUSP  scheme,  Zha  E-CUSP2  scheme, 
which  is  consistent  with  the  characteristic  directions  [1,  3].  The  scheme  has  low 
diffusion  and  is  able  to  capture  crisp  shock  profiles  and  exact  contact  discontinuities. 
The  scheme  is  more  CPU  efficient  since  it  only  uses  the  scalar  dissipation.  In  addition, 
is  is  fairly  straightforward  to  extend  the  new  scheme  to  the  3D  moving  grid  system[27, 
4,  28].  This  is  because  the  grid  velocity  belongs  to  the  convective  terms  in  the  E- 
CUSP  schemes.  The  pressure  term  is  determined  by  the  weighted  average  based  on 
the  wave  eigenvalues  from  downstream  and  upstream.  The  new  E-CUSP  scheme 
is  more  efficient  than  the  Roe  scheme  without  matrix  operation.  For  a  2D  nozzle 
calculation  for  comparison,  the  CPU  time  to  evaluate  the  flux  using  the  new  E-CUSP 
scheme  is  only  about  1/4  of  that  needed  by  the  Roe  scheme  [1], 


3.2  Implicit  Time  Marching  Scheme 

Among  the  researchers  in  the  area  of  3D  time-marching  aeroelastic  analysis  based  on 
Euler/Navier-Stokes  approaches,  Lee-Rausch  and  Batina[29][?]  used  a  three-factor, 
implicit,  upwind-biased  Euler/Navier-Stokes  approach  coupled  with  a  lagged  struc¬ 
ture  solver.  Morton,  Melville  and  Gordnier  et  al.  developed  an  implicit  fully  coupled 
fluid-structure  interaction  model,  which  used  the  Beam-Warming  implicit  approxi¬ 
mate  factorization  scheme  for  the  flow  solver  coupled  with  modal  structural  solver 
[30]  [31]  [14][32].  Liu  et  al.  developed  a  fully  coupled  method  using  Jameson’s  explicit 
scheme  with  multigrid  approach  utilizing  Euler  equations  and  a  modal  structural 
model[13].  Doi  and  Alonso[19]  coupled  an  explicit  Runge-Kutta  multigrid  RANS 
flow  solver  with  a  FEM  structure  solver  to  predict  the  aeroelastic  responses  of  NASA 
Rotor  67  blade. 

In  this  research,  we  have  developed  an  implicit  time  marching  algorithm  using  a 
dual-time  stepping  unfactored  line  Gauss-Seidel  iteration  [5,  2,  4,  28].  The  unfactored 
Gauss-Seidel  iteration  is  unconditionally  stable  and  allows  larger  pseudo  or  physical 
time  steps  than  explicit  method.  It  avoids  the  factorization  error  introduced  by  those 
implicit  approximate  factorization  methods,  such  as  those  used  in  [29]  [30]  [3 1]  [14]  [32] . 
Even  though  the  factorization  error  diminishes  within  each  physical  time  step,  the 
factorization  error  can  limit  the  numerical  stability.  The  linear  stability  analysis 
shows  that  approximate  factorized  method  is  not  stable  for  3D  computation  even 
though  it  is  stable  for  2D  computation. 


3.3  Non-Reflective  Boundary  Conditions 

The  accuracy  of  unsteady  flow  calculations  relies  on  accurate  treatment  of  boundary 
conditions.  Due  to  the  limitation  of  computer  resources,  usually  only  a  finite  com¬ 
putational  domain  is  considered  for  a  flow  calculation.  This  means  that  we  have  to 
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“cut  off”  the  domain  that  is  not  of  our  primary  interest.  However,  the  cut  boundaries 
may  cause  artificial  wave  reflections,  which  may  include  both  physical  waves  and  nu¬ 
merical  waves[33].  Such  waves  may  bounce  back  and  forth  within  the  computational 
domain  and  may  seriously  contaminate  the  solutions  and  produce  misleading  results. 
This  is  particularly  true  for  internal  flows  such  as  the  flows  in  turbomachinery,  in 
which  the  computational  domain  usually  is  confined  very  near  the  solid  objects.  For 
example,  previous  studies  indicated  that  the  different  treatments  of  numerical  per¬ 
turbation  at  upstream  and  downstream  boundaries  can  change  the  compressor  blade 
stall  inception  pattern  [34]  [35]. 

The  currently  often  used  non-reflective  boundary  conditions  for  unsteady  internal 
flows  are  based  on  eigenvalue  analysis  of  linearized  Euler  equations  developed  by 
Giles[36].  However,  Giles’  method  may  only  apply  to  the  inviscid  solutions  which 
require  the  far  field  flow  to  be  uniform  so  that  the  propagation  waves  have  the  Fourier 
mode  shapes.  For  viscous  flows,  the  mean  flow  in  the  downstream  far  field  region  may 
be  non-uniform  due  to  the  airfoil  or  blade  wakes,  which  means  that  there  will  be  no 
Fourier  mode  shapes.  In  addition,  the  inconsistency  of  the  Navier-Stokes  governing 
equations  for  the  inner  domain  and  linearized  Euler  equations  at  far  field  boundary 
may  also  cause  numerical  wave  reflections. 

The  more  rigorous  treatment  of  non-reflective  boundary  conditions  (NRBC)  for 
Navier-Stokes  equations  is  the  one  suggested  by  Poinsot  and  Lele  in  1992  [33]  for  Di¬ 
rect  Numerical  Simulation  of  turbulence.  However,  the  NRBC  given  by  Poinsot  and 
Lele  in  [33]  is  only  for  the  regular  mesh  aliened  with  the  coordinate  axises  in  Carte¬ 
sian  coordinates.  The  explicit  time  marching  scheme  was  used  in  the  calculation  of 
Poinsot  and  Lele.  For  practical  engineering  applications,  the  body  fitted  generalized 
coordinates  are  usually  necessary.  In  2000,  Kim  and  Lee  [37]  made  an  effort  to  extend 
the  NRBC  of  Poinsot  and  Lele  from  the  Cartesian  coordinates  to  generalized  coor¬ 
dinates.  However,  in  their  derivation,  a  flaw  was  made  by  absorbing  the  eigenvector 
matrix  into  the  partial  derivatives,  their  formulations  apply  only  if:  1)  it  is  ID  equa¬ 
tion;  2)  the  eigenvector  matrix  is  constant  in  the  flow  field;  3)  the  partial  differential 
equations  satisfy  Pfaff’s  condition.  For  multidimensional  Navier-Stokes  equations,  all 
these  three  conditions  are  not  satisfied  [38,  39].  Hence,  the  wave  amplitude  vector 
derived  in  [37]  is  erroneous. 

More  recently,  based  on  the  characteristic  approach  of  Poinsot  and  Lele[33],  Bruneau 
and  Creuse  [40]  suggested  a  variation  of  the  approximate  treatment  of  the  incoming 
wave  amplitude  in  the  exit  boundary  conditions  by  assuming  that  the  pressure  and 
velocity  values  will  “convect”  with  time  to  the  location  where  the  phantom  cells  are 
located.  The  results  show  the  method  works  well.  Prosser  and  Schluter  [41]  used 
an  approach  based  on  a  low  Mach  number  asymptotic  expansion  of  the  the  govern¬ 
ing  equations  to  improve  the  specification  of  time  dependent  boundary  conditions. 
With  the  help  of  the  Local  One-Dimensional  Inviscid  (LODI)  relations,  Moureau  et 
al.  [42]  implemented  characteristic  boundary  conditions  for  multi-component  mixtures 
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in  DNS  and  LES  computations  using  a  modified  NRBC  formulation. 

In  this  research,  we  extend  the  NRBC  system  of  Poinsot  and  Lele[33]  from  Carte¬ 
sian  coordinates  to  generalized  coordinates  and  apply  it  numerically  for  unsteady 
calculations  in  an  implicit  time  marching  method.  In  a  finite  difference  or  finite 
volume  approach,  the  governing  equations  are  more  straightforward  to  be  solved  in 
generalized  coordinates,  in  which  a  complex  physical  domain  becomes  a  rectangular 
computational  domain  (for  2-D  case)  or  a  hexahedral  computational  domain  (for  3-D 
case)  with  equal  grid  spacings.  The  moving  grid  effect  can  be  naturally  included 
in  the  generalized  coordinates.  Strictly  speaking,  for  finite  differencing  or  finite  vol¬ 
ume  methods,  only  solving  the  equations  in  generalized  coordinates  can  preserve  the 
accuracy  of  high  order  numerical  schemes. 

In  general,  implicit  methods  permit  a  larger  time  step  and  are  widely  used  for 
many  practical  applications.  To  be  consistent  with  the  implicit  solver  of  the  inner 
domain,  in  this  paper,  the  NRBC  equations  are  implicitly  discretized  and  solved 
simultaneously  in  a  fully  coupled  manner.  Two  numerical  cases  are  tested  in  this 
paper:  a  vortex  propagating  through  a  outflow  boundary  and  a  transonic  inlet-diffuser 
flow  with  shock/boundary  layer  interaction.  The  numerical  results  indicate  that  the 
present  methodology  is  robust  and  accurate. 

The  strategy  is  to  fully  couple  the  flow  and  structural  solver  within  each  time  step 
by  iterating  the  flow  solution  and  structural  deformation.  Since  the  fully  coupled 
fluid-structural  interaction  is  very  CPU  intensive,  this  research  has  developed  high 
efficiency  high  accuracy  CFD  and  structural  algorithms.  The  following  sub-section 
will  outline  the  developed  methodology.  The  detailed  numerical  algorithms  are  de¬ 
scribed  in  next  chapter. 


3.4  Modal  Structural  Solver 

Since  the  full  annulus  of  a  mistuned  rotor  will  be  calculated,  the  nonlinear  3D  Navier- 
Stokes  equations  are  CPU  intensive.  Hence  it  is  very  important  that  the  structural 
solver  is  CPU  efficient  and  accurate.  Based  on  this  consideration,  the  structural  solver 
completely  based  on  the  finite  element  method  for  the  mistuned  bladed  disk  is  not 
favored  due  to  the  high  CPU  cost. 

For  this  proposed  research,  the  structural  solution  for  a  mistuned  bladed  disk  will 
employ  the  SNM  model  in  time  domain(Yang  and  Griffin,  2001)  [43]  developed  under 
the  support  of  the  GUIde  Consortium  (Government,  Universities,  and  Industry).  The 
SNM  model  uses  a  subset  of  tuned  bladed  disk  modes  to  represent  the  vibration  of  a 
mistuned  bladed  disk.  It  is  verified  that  the  SNM  model  is  both  numerically  accurate 
and  computationally  efficient  (Srinivasan,  1999)  [44]. 

Other  well-recognized  models  for  mistuned  bladed  disks  include  TURBO  REDUCE 
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(Kruse  and  Pierre,  1996) [45]  and  MISTRESS  (Petrov  et.  al.,  2000)  [46].  TURBO 
REDUCE  utilizes  a  Component  Mode  Synthesis  (CMS)  technique  that  accounts  for 
blade  frequency  mistuning  only.  MISTRESS  and  SNM,  on  the  other  hand,  employ 
the  Modal  Reduction  (MR)  technique  which  allows  general  mistuning  of  the  structure 
including  mass  and  stiffness  variations.  Since  it  is  our  desire  to  develop  a  methodol¬ 
ogy  that  can  be  applied  to  general  mistuning  problems,  Modal  Reduction  technique 
(MISTRESS  and  SNM)  became  the  method  of  our  choice.  A  detailed  comparison 
of  CMS  and  MR  techniques  for  mistuned  bladed  disk  vibration  is  documented  by 
Moyroud  et  al,  2002  [47]. 

We  selected  SNM  over  MISTRESS  based  on  the  following  three  reasons: 

1.  The  fluid-structural  interaction  problem  requires  the  calculation  of  the  vibration 
of  all  nodes  on  the  surfaces  of  the  airfoils  of  a  bladed  disk. 

2.  MISTREE  uses  the  receptance  method  which  is  only  efficient  when  the  number 
of  nodal  solutions  calculated  is  limited  (e.g.,  the  study  of  the  vibration  of  a  mistuned 
bladed  disk  with  friction  dampers  only  interests  in  the  vibration  at  a  limited  number  of 
friction  joints.)  Using  MISTRESS  for  our  study  would  be  computationally  expensive 
since  the  number  of  nodes  on  the  airfoil  surfaces  of  a  bladed  disk  can  range  from 
10,000  to  100,000. 

3.  SNM  uses  a  limited  set  of  tuned  modes  to  represent  the  vibration  of  a  mistuned 
bladed  disk.  Its  efficiency  solely  depends  on  the  number  of  tuned  modes  of  choice. 
The  typical  number  of  modes  needed  to  represent  an  industrial  bladed  disk  is  in  the 
order  of  100  (Srinivasan,  1999)  [44]. 

Considering  CFD  calculation  is  very  CPU  intensive,  using  SNM  to  simulate  the 
response  of  mistuned  bladed  disks  becomes  particularly  appealing  and  essential  to 
make  the  simulation  of  the  fully  coupled  fluid-structural  problem  possible. 


3.5  High  Performance  Computing 

The  parallel  computing  capability  based  on  SPMD  (Single  Program  Miltiple  Data) 
is  implemented  in  our  code  to  reduce  the  wall  clock  calculation  time.  The  reduction 
of  wall  clock  time  by  parallel  computing  is  essential  and  necessary. 
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Chapter  4 

Discretization  Schemes 

4.1  Flow  Governing  Equations 


The  governing  equations  for  the  flow  field  computation  are  the  Reynolds-Averaged 
Navier-Stokes  equations  (RANS)  with  Favre  mass  average  which  can  be  transformed 
to  the  generalized  coordinates  and  expressed  as: 

dQ'  dE'  OF'  dGr  1  (OE'  8F'V  dG'\ 

■y+sr+a?+^  =  s(^f+  (41) 

where  Re  is  the  Reynolds  number  and 


Q'  =  j  (4.2) 

E'  =  J&Q  4-  &E  +  £yF  +  £2G)  =  j(6Q  +  E")  (4.3) 

F'  =  j(Vt  Q  +  rjxE  +  VyF  +  VzG)  =  j(rjt  Q  +  F")  (4.4) 

G'  =  j((tQ  +  CxE  +  CyF  +  CG)  =  i(CtQ  +  G")  (4.5) 

Ey  =  -j(£xEv  +  £yFv  +  £2Gv)  (4.6) 

F(,  =  —  (r)xEv  +  T7yFv  +  7y2Gv)  (4.7) 
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G'v  =  ~(CxEv  +  CyFv  +  C2Gv)  (4.8) 

where  the  variable  vector  Q,  and  inviscid  flux  vectors  E,  F,  and  G  are 
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The  E",  F",  and  G"  are  the  inviscid  fluxes  at  the  stationary  grid  system  and  are: 


Eff  =  £BE  +  &F  +  &G, 
F"  =  ??XE  +  rjyF  +  r)zG, 


G"  =  CxE  +  CyF  +  C,G, 


and  the  viscous  flux  vectors  are  given  by 


Ev  = 
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Txx  ~  (m”u" 

f  xy  ~  PU"V" 


Txz  —  pu"w" 


Qz 


,FV  = 


Tyx  -  pv"u" 


Tyy  -  PV"V" 


Tyz  —  PV  W 
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,  Gv  = 


fzx  —  pw"u" 
fzy  —  pw"v" 


Tzz  —  pw"w" 


Qz 


In  above  equations,  p  is  the  density,  u,  v,  and  w  are  the  Cartesian  velocity  compo¬ 
nents  in  x,  y  and  z  directions,  p  is  the  static  pressure,  and  e  is  the  total  energy  per  unit 
mass.  The  overbar  denotes  the  Reynolds-averaged  quantity,  tilde  and  double-prime 
denote  the  Favre  mean  and  Favre  fluctuating  part  of  the  turbulent  motion  respec¬ 
tively.  All  the  flow  variable  in  above  equations  are  non-dimensionlized  by  using  the 
freestream  quantities  and  a  reference  length  L. 

Let  subscript  1,  2  and  3  represent  the  coordinates,  x,  y,  and  z,  and  use  Einstein 
summation  convention,  the  shear-stress  and  Qx,  Qy,  Qz  terms  in  non-dimensional 
forms  can  be  expressed  in  tensor  form  as 


2  _  diik  r  _ ,  dui  duj . 


(4.9) 


Qi  =  Uj{Tij  -  pu"u")  -  ( Qi  +  CvpT"u") 


(4.10) 
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where  the  mean  molecular  heat  flux  is 


_  p  da 2 

^  (7  —  l)Pr  dxi 

The  molecular  viscosity  p  =  p(T)  is  determined  by  Sutherland  law,  and  a  = 
ClRToo  is  the  speed  of  sound.  The  equation  of  state  closes  the  system, 

pe  =  7  — — r  +  ~p(u7  +  V2  +  W2)  +  k  (4.12) 

(7  -  1J  2 

where  7  is  the  ratio  of  specific  heats,  k  is  the  Favre  mass-averaged  turbulence  kinetic 
energy.  The  turbulent  shear  stresses  and  heat  flux  appeared  in  above  equations  are 
calculated  by  Baldwin-Lomax  model[48].  The  viscosity  is  composed  of  p  +  pt,  where 
/i  is  the  molecular  viscosity  and  pt  is  the  turbulent  viscosity  determined  by  Baldwin 
Lomax  model.  For  a  laminar  flow,  the  pt  is  set  to  be  zero. 


4.2  Time  Marching  Scheme 

The  time  dependent  governing  equation  (4.1)  is  solved  using  the  control  volume 
method  with  the  concept  of  dual  time  stepping  suggested  by  Jameson  [49].  A  pseudo 
temporal  term  ^  is  added  to  the  governing  equation  (4.1).  This  term  vanishes  at 
the  end  of  each  physical  time  step,  and  has  no  influence  on  the  accuracy  of  the  solu¬ 
tion.  However,  instead  of  using  the  explicit  scheme  as  in  [49],  an  implicit  pseudo  time 
marching  scheme  using  line  Gauss-Seidel  iteration  is  employed  to  achieve  high  CPU 
efficiency.  For  unsteady  time  accurate  computations,  the  temporal  term  is  discretized 
implicitly  using  a  three  point,  backward  differencing  as  the  following 

<9Q  _  3Qn+1  -  4Q"  +  Qn_1 

dt  2Af  ^  ' 

Where  n  is  the  time  level  index.  The  pseudo  temporal  term  is  discretized  with 
first  order  Euler  scheme.  Let  m  stand  for  the  iteration  index  within  a  physical  time 
step,  the  semi-discretized  governing  equation  (4.1)  can  be  expressed  as 


l(i  +  = 


3Qn+l,m  _  4Q n  +  Qn-1 

2A t 


(4.14) 


where  the  At  is  the  pseudo  time  step,  R  is  the  net  flux  going  through  the  control 
volume, 
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R  “  ~vfP'  ~  +  &  -  +  <G'  -  (415> 

where  V  is  the  volume  of  the  control  volume,  s  is  the  control  volume  surface  area  vec¬ 
tor.  Equation  (4.14)  is  solved  using  the  unfactored  line  Gauss-Seidel  iteration.  Two 
line  sweeps  in  each  pseudo  time  steps  are  used,  one  sweeps  forward  and  the  other 
sweeps  backward.  The  alternative  sweep  directions  are  beneficial  to  the  information 
propagation  to  reach  high  convergence  rate.  Within  each  physical  time  step,  the  so¬ 
lution  marches  in  pseudo  time  until  converged.  The  method  is  unconditionally  stable 
and  can  reach  very  large  pseudo  time  step  since  no  factorization  error  is  introduced. 


4.3  The  Zha  E-CUSP  Scheme[l,  2] 


To  clearly  describe  the  formulations,  the  vectors  Q  and  E'  in  Eq.  (4.1)  are  given 
below: 
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(4.16) 


U  is  the  contravariant  velocity  in  f  direction  and  is  defined  as  the  following: 


U  =  &  +  fxtt  +  £yv  +  £zw 


U  is  defined  as: 


(4.17) 


U  =  U-&  (4.18) 

The  Jacobian  matrix  A  is  defined  as: 

<9E  -  „  „ 

A  =  ~  =  TAT"1  (4.19) 

where  T  is  the  right  eigenvector  matrix  of  A,  and  A  is  the  eigenvalue  matrix  of  A  on 
the  moving  grid  system  with  the  eigenvalues  of: 


(U  +  C,U-C,  U,  U,  U) 


(4.20) 
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where  C  is  the  speed  of  sound  corresponding  to  the  contravariant  velocity: 

C  =  Cy/g  +  %  +  g  (4.21) 

and  where  c  =  y/jRT  is  the  physical  speed  of  sound. 

Due  to  the  homogeneous  relationship  between  Q  and  E,  the  following  formulation 
applies: 


E  =  AQ  =  TAT_1Q  (4.22) 

In  an  E-CUSP  scheme,  the  eigenvalue  matrix  is  split  as  the  following: 
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The  grid  velocity  term  £t[I]  due  to  the  moving  mesh  is  naturally  included  in  the 
convective  term,  U,  as  given  in  Eq.  (4.43).  Therefore,  Eq.  (4.22)  becomes: 
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(4.24) 


where  Ec  and  Ep  are  namely  the  convective  and  pressure  fluxes.  As  shown  above, 
the  way  of  splitting  the  total  flux  into  convective  and  pressure  fluxes  in  an  E-CUSP 
scheme  is  purely  based  on  the  analysis  of  characteristics  of  the  system.  As  shown  in 
Eq.  (4.24),  the  convective  flux  has  the  upwind  characteristic  U  and  is  only  associated 
with  the  convective  velocity.  The  pressure  flux  has  a  downwind  and  an  upwind 
characteristic  and  it  completely  depends  on  the  propagation  of  an  acoustic  wave. 
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The  Zha  E-CUSP2  scheme  is  based  on  the  E-CUSP  scheme  suggested  by  Zha  and 
Hu  [1],  which  is  extended  to  a  moving  mesh  system  by  the  following: 
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T. 

U  P(U-Cx)J 

(4.25) 


where 


(pU)h  =  (pLU++pRUn) 


(  1  \ 
u 
v 
w 

\  i  ) 


Cl  =  \(CL  +  CR) 


UL  u„ 

2  2 


.  Mi  +  IM/.I 


0?  =  Ci{—  V"""  +  +  1)!  -  V— ')? 


Ml  +  |Ml| 


(4.26) 

(4.27) 

(4.28) 

(4.29) 

(4.30) 


.Mr  —  \Mr\ 


OIL 


2 

2  (p/p)l 


0-R  =  +  aR[-\(MR  -  If  - 


Mr  —  \Mr\ 


( P/p)l  +  {P/p)r 


o:r 


2  (P/p)r 


( P/p)l  +  (p/p)r 


P±  =  i(M±l)2(2TM)±aM(M2-  l)2,  a  =  ^ 


(4.31) 

(4.32) 

(4.33) 


C  =  C-£t 


(4.34) 
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C\  =  \(Cl  +  Cr)  (4.35) 

Please  note  that,  in  the  energy  equation  of  the  pressure  splitting,  U  and  C  are 
used  instead  of  U  and  C.  The  term  C  is  constructed  by  taking  into  account  the  effect 
of  the  grid  velocity  so  that  the  flux  will  transit  from  subsonic  to  supersonic  smoothly. 
When  C  =  0,  Eq.  (4.25)  naturally  returns  to  the  one  for  a  stationary  grid. 

For  supersonic  flow,  when  Ul  >  C,  Ei  =  E/,;  when  Ur  <  —C.  Ei  =  E/j. 

2  2 

4.3.1  Numerical  Dissipation 

The  low  numerical  dissipation  at  stagnation  is  important  to  accurately  resolve  wall 
boundary  layers.  An  upwind  scheme  can  be  written  as  a  central  differencing  plus  a 
numerical  dissipation. 

To  analyze  the  numerical  dissipation  at  stagnation,  the  ID  Euler  equation  is  used 
as  the  example.  Assuming  u  =  0,  the  numerical  dissipation  vector  of  the  new  E-CUSP 
scheme  at  stagnation  is: 


D  =  — -f 


a  i 
2 


/  0  \ 
0 

\  6P 


where 


6p  =  pR-  pL 


The  numerical  dissipation  of  the  Roe  scheme  at  stagnation  is: 


(4.36) 


(4.37) 


a± 

n  ^  _ _ 2 

uRoe-  2(7_1) 


( (7 

V 


l)/a\5p  \ 
2 

0 

dp 


(4.38) 


where  the  "stands  for  the  Roe’s  average[50]. 

Comparing  eq.(4.36)  and  (4.38),  it  can  be  seen  that  the  numerical  dissipation  of 
the  new  E-CUSP  scheme  for  the  continuity  equation  vanishes  at  u  =  0  while  the  Roe 
scheme  has  the  non-vanishing  dissipation.  For  the  energy  equation,  the  two  schemes 
have  equivalent  dissipation.  For  ideal  gas  with  the  7  =  1.4,  the  coefficient  of  the 
Roe  scheme  energy  dissipation  term  is  2.5  times  larger  than  that  of  the  new  E-CUSP 
scheme. 


In  conclusion,  even  though  there  is  one  non-vanishing  numerical  dissipation  term 
in  the  energy  equation  for  the  new  E-CUSP  scheme,  the  overall  numerical  dissipation 
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of  the  new  E-CUSP  scheme  is  not  greater  than  that  of  the  Roe  scheme.  The  Roe 
scheme  is  proved  to  be  accurate  to  resolve  wall  boundary  layers[51].  It  is  hence 
expected  that  the  new  E-CUSP  scheme  should  also  have  sufficiently  low  dissipation 
to  accurately  resolve  wall  boundary  layers.  This  is  indeed  the  case  shown  by  the 
numerical  experiment  for  a  flat  plate  boundary  layer. 


4.3.2  Zha  E-CUSP2  Scheme[3] 

The  original  Zha-Hu  scheme  is  proved  to  have  low  diffusion  and  is  able  to  capture 
crisp  shock  wave  profile  and  exact  contact  surface[l].  However,  the  scheme  is  found 
to  have  temperature  oscillations  near  the  solid  wall  region  when  the  grids  is  skewed. 
Therefore,  the  scheme  used  in  the  present  study  is  the  modified  version  scheme,  Zha 
CUSP2  scheme[52].  In  this  scheme,  the  total  enthalpy  instead  of  the  static  pressure 
is  used  to  calculate  the  numerical  dissipation  coefficients  for  the  energy  equation  as 
below: 


OiL 


2  (H/p)l 


(H/p)L  +  {H/p)R' 


aR 


2( H/p )j 


(H/p)L  +  (H/p)R 


(4.39) 


Note  that  Equation  (4.39)  is  only  used  for  the  energy  equation.  For  the  continuity 
and  momentum  equations,  Equation  (4.32)  is  still  used  as  the  smoothing  coefficient. 


4.4 


Roe’s  Riemann  Solver  on  Moving  Grid  System[4, 
5] 


The  Roe’s  Riemann  solver  is  also  implemented  in  the  solver  as  a  benchmark  scheme 
to  compare  the  results.  Roe  scheme  is  recognized  as  having  very  low  diffusion  and 
can  capture  exact  shock  and  contact  discontinuities.  In  present  study,  the  original 
Roe  scheme  is  extended  to  moving  grid  system  as  the  following,  for  example,  in  £ 
direction: 


E'+|  =  -j[ E"(Ql)  +  E"(Qr)  +  QL£tL  +  Qr6*  -  |A|(Qr  -  QL)]i+i  (4.40) 

where  Ql  and  Qr  are  the  reconstructed  variables  to  the  left  and  right  sides  of  the 
cell  face,  a°d  Cr  are  the  reconstructed  grid  velocity  component  in  £  direction  to 
the  left  and  right  sides  of  the  cell  interface  i  +  A  is  the  Jacobian  matrix,  A  =  Of— 
and  it  takes  the  form  as  A  =  TAT-1,  T  is  the  right  eigenvector  matrix  of  A,  A  is 
the  eigenvalue  matrix  of  A,  and 
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A  =  TAT"1  (4.41) 

where  A  is  the  eigenvalue  matrix  on  moving  grid  system  with  the  eigenvalues  of 

(U  +  C,U-C,U,U,  U)  (4.42) 

where  U  is  the  contravariant  velocity  in  £  direction  on  moving  grid, 

u  =  £t  +  txU  +  €yV  +  £ZU)  (4.43) 

C  is  the  speed  of  sound  corresponding  to  the  contravariant  velocity: 

c  =  Cy/g  +  %  +  Q  (4.44) 

where  c  =  y/j. RT  is  the  physical  speed  of  sound.  The  ~  stands  for  the  Roe- averaged 
quantities.  For  example, 

it  =  (&l  +  &r\J Pr/ Pl)  /  (1  +  \Jpr/pl)  (4.45) 

It  can  be  proved  that  the  eigenvector  matrix  T  has  exactly  the  same  form  as 
the  one  without  moving  grid.  The  only  difference  between  the  moving  grid  and  the 
stationary  grid  system  is  that,  for  the  moving  grid  system,  the  contravariant  velocity 
in  the  eigenvalues  contains  the  grid  velocity  as  given  in  Equation  (4.43).  It  is  hence 
straightforward  to  extend  the  code  from  a  stationary  grid  system  to  the  moving  grid 
system  using  Roe  scheme  without  major  change. 

The  grid  velocity  is  evaluated  at  the  center  of  each  cell  and  is  determined  by 
the  averaged  value  that  counts  the  movement  of  the  eight  vertexes  if  hexahedral 
control  volumes  are  used.  The  grid  velocity  is  reconstructed  with  3rd  order  MUSCL 
differencing. 


4.5  Conventional  Boundary  Conditions 

Two  sets  of  boundary  conditions  are  developed  in  this  research.  The  first  set  is  the 
conventional  boundary  conditions.  The  second  set  is  the  non-reflective  boundary 
conditions  to  be  described  in  Chapter  5. 

The  conventional  boundary  conditions  used  for  both  the  steady  state  and  unsteady 
calculation  are  as  follows: 

(1)  Inlet  boundary  conditions:  The  far  field  boundary  is  divided  into  inlet  and 
outlet  boundaries.  On  inlet  boundary,  it  is  assumed  that  the  streamwise  velocity  u 
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is  uniform,  transverse  velocity  v  =  0,  and  spanwise  velocity  w  =  0.  Other  primitive 
variables  are  specified  according  to  the  freestream  condition  except  the  pressure  which 
is  extrapolated  from  interior. 

(2)  Outlet  boundary  conditions:  All  the  flow  quantities  are  extrapolated  from 
interior  except  the  pressure  which  is  set  to  be  its  freestream  value. 

(3)  Solid  wall  boundary  conditions:  At  moving  boundary  surface,  the  no-slip  con¬ 
dition  is  enforced  by  extrapolating  the  velocity  between  the  phantom  and  interior 
cells, 


uo  =  2xb-uu  v0  =  2yb-v1,  w0  =  2zb-wi  (4.46) 

where  Uo,  v0  and  w0  denote  the  velocity  at  phantom  cell,  u\,  v\  and  W\  denote  the 
velocity  at  the  1st  interior  cell  close  to  the  boundary,  and  ub,  vb  and  wb  are  the  velocity 
on  the  moving  boundary. 

If  the  wall  surface  is  in  rj  direction,  the  other  two  conditions  to  be  imposed  on 
the  solid  wall  are  the  adiabatic  wall  condition  and  the  inviscid  normal  momentum 
equation  [30]  as  follows, 

?  =  °,  f?  =  -  (  ~2~,  ^2-  --2  )  (VxXb  +  vyyb  +  rjzzb)  (4.47) 

dr]  dy  \Vi  +  Vy  +  mJ 

4.6  Moving/Deforming  Grid  Systems 

In  the  fully-coupled  computation,  the  remeshing  is  performed  in  each  iteration.  There¬ 
fore,  a  CPU  time  efficient  algebraic  grid  deformation  method  is  employed  in  the  com¬ 
putation  instead  of  the  commonly-used  grid  generation  method  in  which  the  Poisson 
equation  is  solved  for  grid  points.  For  clarity,  the  remeshing  procedure  for  2D  cases  is 
sketched  in  Figure  12.131.  This  grid  deformation  procedure  is  designed  in  such  a  way 
that  the  far-field  boundary  (j=jlp)  is  held  fixed,  and  the  grids  on  the  wing  surface 
(j=l)  moves  and  deforms  following  the  instantaneous  motion  of  the  wing  structure. 
After  the  new  wing  surface  is  determined,  two  components  of  the  displacement  vector 
at  wing  surface  node  dx i  and  dy\  can  be  calculated  accordingly.  First,  the  length  of 
each  segment  along  the  old  mesh  line  is  estimated  as: 

Sj  =  i  +  y/{xj  -  Xj-x)2  +  fa  -  yj- 1)2  (j  =  2,  •  •  • ,  jlp)  (4.48) 

where  S\  =  0  and  the  displacement  vectors  at  wing  surface  node  (  dx\,  dy\  )  and  at 
the  far-field  boundary  (  dxjip,  dy}lp  )  are  known.  Then  the  grid  node  points  between 
the  wing  surface  and  the  far-field  boundary  can  be  obtained  by  using  following  linear 
interpolation: 
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,  dxjip  dx\  dxiSjip  dxjipSi 

CtX j  —  S  j  I 

Sjlp  Ip  ®1 

j. .  dUjip  dy\  _  dtfoSjip  dyjipsi 

(*Vj  —  T 

Sjlp  Si  Sjlp  S  1 


(4.49) 

(4.50) 


This  simple  remeshing  strategy  is  proved  to  be  robust  for  all  the  cases  investigated 
in  present  study.  By  monitoring  the  accuracy  criterion  y+ ,  it  is  shown  that  the  method 
can  maintain  the  initial  grid  quality  and  keep  almost  the  same  mesh  distribution 
around  the  wing  surface. 

For  3D  case,  the  Equation  (4.51)  becomes 


sj  =  Sj—i  +  yj (xj  ~  Xj-i)2  +  {yj  ~  t/j-i)2  +  [zj  -  Zj- 1)2  0  =  2,---,  jlp)  (4.51) 

and  one  more  equation  is  added  to  determined  the  z  component  of  displacement 
vector: 


dzj  = 


dzjip  dz\  dz\  s  jip  dzjipSi 

gj  _j_ 

$jlp  Sjlp 


(4.52) 


4.7  Geometric  Conservation  Law 


It  was  pointed  out  by  Thomas  et  al.  [53]  that  due  to  the  mixed  temporal  and  spa¬ 
tial  derivatives  after  discretization,  an  additional  term  appears,  which  theoretically 
equals  to  zero  but  numerically  still  remains.  Consequently  numerical  error  could  be 
introduced  in  the  discretized  form  of  the  equations  of  the  flow  motion  if  this  term 
is  neglected.  In  order  to  reduce  or  avoid  this  error,  the  geometric  conservation  law 
needs  to  be  enforced.  In  other  words,  the  following  additional  term  should  be  added 
to  the  right-hand  side  of  the  equations  as  a  source  term: 


S  =  Q 


(4.53) 


To  implement  this  option  in  the  flow  solver,  the  source  term  is  then  linearized  such 
that 


C"+l 


=  s.,  +  |LAcr. 


(4.54) 


As  has  been  observed  in  ref.  [5],  the  overall  performance  of  this  numerical  supple¬ 
ment  is  beneficial  with  very  little  CPU  time  cost. 
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Chapter  5 

Non-Reflective  Boundary 
Conditions 


5.1  Characteristic  Form  of  the  Navier-Stokes  Equations[6] 


The  characteristic  form  of  the  Navier-Stokes  equations  in  the  generalized  coordinates 
will  be  solved  to  determine  the  non-reflective  boundary  conditions  at  the  phantom 
cells.  To  describe  the  derivation  process,  the  £  direction  will  be  taken  as  an  exam¬ 
ple.  For  the  other  two  directions,  the  formulations  can  follow  the  same  procedure 
and  the  general  formulations  are  given  in  the  appendix.  Based  on  the  strategy  of 
Thompson[39]  and  Poinsot  and  Lele[33],  the  Navier-Stokes  equations  are  expressed 
first  using  primitive  variables  as  the  following: 

M^  +  A-M^  +  B-M^  +  C-M^  =  RV  (5.1) 

at  d£  or]  d£ 

where  A,B,  C  are  the  Jacobian  matrix 

dR  dF  dCT 

dQ'1  dQ'’  dQ' 

where  Rv  is  the  viscous  vector  on  the  right  hand  side  of  the  Navier-Stokes  equations, 
(  Equation  (4.1)  ),  q  is  the  primitive  variable  vector: 


(  P\ 


\  P  ) 


(5.3) 
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M  is  the  Jacobian  matrix  between  the  conservative  variables  and  primitive  variables 

/I  0  0  0  0  \ 

an;  up  0  0  0 

M  =  -^-  =  v  0  p  0  0  (5.4) 

dq  w  0  0  p  0 

pu  pv  pw  ^  J 

where  $  =  -^(u2  +  v2  +  w2). 

Equation  (5.1)  can  be  further  expressed  as: 


Where 


dq  dq  ,  dq  dq 

¥  +  a'5f+bs^  +  cac=  ^ 


(5.5) 


a  =  M_1AM,  b  =  c  =  M_1CM  (5.6) 


0  0  0  0  \ 

j  0  0  0 

o  }  0  0  (5.7) 

0  0  }  0 

u(y  —  1)  —  v(y  —  1)  —w(  7  —  1)  7  —  1 

Matrix  a,  b,  c  have  the  same  eigenvalues  as  Jacobian  matrix  A,  B,  C.  In  £  direction, 

'  U  pC  p(y  pC  0  \ 

0  U  0  0  ^ 

a  =  0  0  U  0  “  (5.8) 

0  0  0  U 

p 

\  0  7 pC  7 pC  t pC  u 

■where  U  =  Cu  +  Cv  +  Cw-  Matrix  a  can  also  be  expressed  as 

a  =  PAP"1  (5.9) 

where  A  is  the  eigenvalue  matrix,  P  is  eigenvector  matrix  of  a,  and  P_1  is  the  inverse 
of  P.  They  are  given  as  the  following 
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U  0  0  0  0 

0  U  0  0  0 

0  0  U  0  0 

000  U+C  0 
0  0  0  0  u-c 


(5.10) 


£ y  ^ z  ^  OL 

0  -i,  iy  UV2  -1/V2 

1  o  -fx  iy/y/2  -iy/y/2 

-iy  L  0  l/y/2  -UV2 

0  0  0  ac2  ac2 


(5.11) 


Zx  o 

iy  -L 


-iy  -i*/<? 

L  -iyt<? 

o  -6/c2 


o  (X/V2  iy/s/2  l/y/2  0 

o  -L/V2  -iy/y/2  -l/y/2  (3 


(5.12) 


where  C  =  c|  V£|,  |V£|  =  fa2  +  tv2  +  £*,  L  =  &/{?&  |y  =  &  =  &/|V£|, 

a  =  p/y/2c,  0  =  \jy/2pc ,  c  is  the  speed  of  sound  and  determined  by  c  =  y/'yRT. 

The  Navier-Stokes  equation,  Equation  (5.5)  then  can  be  expressed  as: 


£  +  PAP--«3  +  b£+«S_M->I«. 

o£  dp  d( 


(5.13) 


+ AP_I| + p~'b| +  p'Ic|  =  p‘lM_lR*  <5-14> 

This  is  the  characteristic  form  of  the  Navier-Stokes  equations  in  /  direction.  Define 
vector  C  as: 


— -i 

The  Navier-Stokes  equations  (Equation  (5.14))  are  then  expressed  as: 


(5.15) 


P -’(j!  +  £  +  P“'b^3  +  p-‘c25  =  P-1M_1RV 

at  on  dC 


(5.16) 


Vector  C  is  given  as  the  following: 
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C  — 


(  Cl  \ 

/ 

c3 

= 

£4 

Us  / 

V 

"MW + MW  -  MW 

lUl)  +  14(7)  c^v, 
^4(5)+4|(j)-i#(7)-&|v, 


%£«)] 


(5)1 


+ c)[4|0)  +  Mil)  +  Mil) + /*4(5)1 


(5.17) 


The  vector  £  is  the  amplitude  of  the  characteristic  waves.  If  assume  (x  =  l,£y  = 
£z  =  0,  Equation  (5.17)  returns  to  the  corresponding  formulations  in  x-direction  of 
the  Cartesian  coordinates. 

As  pointed  out  in  [38,  39],  for  multi-dimensional  Navier-Stokes  flow  equations, 
Equation  (5.14),  the  matrix  P-1  can  not  be  absorbed  into  the  partial  derivatives 
because  the  flow  equations  does  not  satisfy  Pfaff’s  condition  and  the  matrix  can  not 
be  treated  as  constants.  In  other  words,  it  is  incorrect  to  express  the  characteristic 
form  of  the  Navier-Stokes  equations  in  the  form  given  in  [37]  (page  2042)  as: 


dR  dR  j 

_  +  A_+  =  P>s; 


(5.18) 


The  local  one-dimensional  wave  amplitude  defined  in  [37]  following  Equation  (5.18) 
is  therefore  also  erroneous. 

To  be  consistent  with  the  governing  equations  of  the  flow  field  within  inner  domain 
and  facilitate  programming,  it  is  desirable  to  express  Equation  (5.16)  in  terms  of 
conservative  variables.  Multiply  Equation  (5.16)  by  matrix  M  •  P,  the  characteristic 
Navier-Stokes  equations  expressed  in  terms  of  conservative  variables  in  £  direction  is: 


dQ' 

dt 


<9F' 

+  MPT  +  —  + 
or) 


dG' 

sc 


1  (dE'v  d¥'v  dG'A 
Re\  d£  dr)  dC  ) 


Define  vector  d  as 


(5.19) 


(  dj  ) 

^  £xTi  +  £j/£ 2  T  3  T  ®(£ 4  T  £5)  ^ 

dz 

~£zPt  +  +  ^(£4  —  £5) 

a 

II 

TJ 

h> 

II 

d3 

= 

£z£  1  -  £x^3  +  ^(£4  —  £5) 

d4 

\  ds  ) 

-£y£l  -  1x^2  +  ^(£4  -  £5) 

<  ac2(£4  +  £5)  / 

Define  vector  V  as: 
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V  =  Md  = 


(  d>  \ 

udj  +  pdg 

vd2  +  pd3 

wdi  +  pd4 

K  \{u2  +  v2  +  w2)d i  +  pudz  +  pvd3  +  pwd4  +  d5  j 


(5.21) 


Finally  the  characteristic  form  Navier-Stokes  equations  in  conservative  form  and 
generalized  coordinates  in  £  direction  can  be  expressed  as: 

m  v  ,  9F'  0Gf  1  fdE'  dK  ,dG\r\ 

6t  dr)  d(  Re\  d£  dr)  d(  ) 

Equation  (5.22)  will  be  solved  to  determine  the  non-reflective  boundary  conditions 
in  £  direction.  The  Navier-Stokes  equations  in  generalized  coordinates  and  their 
characteristic  forms  in  rj  and  (  directions  can  be  obtained  straightforwardly  following 
the  symmetric  rule  and  are  given  in  the  appendix. 

By  neglecting  the  transverse  and  viscous  terms  in  Equation  (5.22),  the  Local  One- 
Dimensional  Inviscid  (LODI)  relation  [33]  in  generalized  coordinates  is 

^  +  £>  =  0  (5.23) 

The  LODI  relation  may  be  used  to  estimate  the  amplitudes  of  the  characteristic 
waves  at  boundaries.  Numerical  results  show  that  the  LODI  relations  works  well  for 
the  boundaries  where  the  flow  fields  are  smooth  or  uniform,  and  hence  the  transverse 
and  viscous  terms  are  small  or  negligible.  For  those  boundaries  where  the  transverse 
and  viscous  terms  are  significant,  the  LODI  relations  may  perform  poorly. 


5.2  Non-Reflective  Boundary  Conditions [6] 

Following  the  strategy  suggested  by  Poinsot  and  Lele[33],  the  characteristic  boundary 
conditions  for  Navier-Stokes  equations  can  be  implemented  based  on  Equation  (5.22). 
In  the  present  study,  Equation  (5.22)  is  solved  implicitly  at  the  phantom  cells  in  a 
fully  coupled  manner  with  the  Navier-Stokes  equations  governing  the  inner  flow  field. 
For  unsteady  solutions,  the  dual  time  stepping  method  is  used.  The  semi-discretized 
equation  for  Equation  (5.22)  is: 

V  1  1.5\  j  fdRbc\n+1’m  (&D\ 

( Ar  ^  At)  1  (  OQJ  +\dQj 


6Qn+1’m+ 1 
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nn+l,m 

nbc 


-  T>n+l 


,771 


3Qn+l,m 


-  4  or  +  Qn~l 
2A  t 


(5.24) 


where 


R-  -  +  <F'  -  +  <G'  ~  iSG”)k]'ds  <5-25) 

Compare  Equations  (5.25)  and  (4.15),  it  is  noticed  that  in  Rbc,  there  is  no  E'  flux, 
which  is  replaced  by  vector  V.  V  is  treated  as  a  source  term. 

Before  proceeding  to  the  further  analysis,  some  notations  need  to  be  defined.  For 
the  finite  volume  method  used  in  the  present  study,  a  row  of  phantom  cells  are  used 
outside  of  the  boundary.  The  boundary  conditions  are  enforced  by  assigning  values 
to  the  primitive  variables  at  those  phantom  cells.  All  the  variables  marked  by  the 
subscript  ‘o’  are  for  phantom  cells.  The  variables  at  the  interior  cells  adjacent  to  a 
boundary  are  denoted  by  subscript  ‘i\ 

Equation  (5.22)  provides  the  set  of  governing  equations  for  NRBC,  but  the  way  to 
implement  the  NRBC  is  not  unique.  The  following  is  the  method  used  in  this  study 
and  should  not  be  considered  as  the  only  feasible  method. 


5.2.1  Supersonic  outflow  boundary  conditions 

For  supersonic  flow  at  exit,  all  the  eigenvalues  in  Equation  (5.10)  are  positive  and  the 
disturbance  propagates  from  inner  domain  to  outside.  The  wave  amplitude  vector, 
Equation  (5.17)  is  evaluated  using  one  side  upwind  differencing.  For  supersonic  flow 
at  exit,  using  simple  extrapolation  may  not  generate  physical  wave  reflection,  but  may 
still  generate  numerical  wave  reflection [33].  Solving  Equation  (5.24)  would  achieve 
a  more  accurate  non-reflective  boundary  conditions  for  the  supersonic  flow.  For  su¬ 
personic  flow,  the  exit  boundary  conditions,  pQ,  pua,  pva,  pwa  and  peQ  are  completely 
determined  by  solving  the  Navier-Stokes  equations  in  the  characteristic  form. 

To  evaluate  the  derivatives  in  vector  C,  either  the  first  order  or  second  order  upwind 
differencing  may  be  used.  For  the  present  study,  all  the  partial  derivatives  in  vector 
£  are  calculated  by  first  order  upwind  differencing. 

5.2.2  Subsonic  outflow  boundary  conditions 

For  subsonic  flow  at  exit,  the  eigenvalue  U  —  C  is  negative  and  the  disturbance 
propagates  into  the  domain  from  outside.  £x  to  £4  can  be  still  calculated  by  one-side 
upwind  differencing.  However,  £5  corresponding  to  the  eigenvalue  of  U  —  C  must  be 
treated  differently.  The  commonly-used  method  to  provide  a  well  posed  boundary 
condition  is  to  impose  p  =  p^  at  the  outflow  boundary.  This  treatment  however 
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will  create  acoustic  wave  reflections,  which  may  be  diffused  and  eventually  disappear 
when  the  solution  is  converged  to  a  steady  state  solution.  For  unsteady  flows,  the 
wave  reflection  may  contaminate  the  flow  solutions.  To  avoid  wave  reflections,  the 
following  soft  boundary  condition  was  suggested  by  Rudy-Strikwerda[54]  and  used  by 
Poinsot-Lele[33]. 


C5  =  IC(p-pe)  (5.26) 

where  1C  is  a  constant  and  is  determined  by  1C  =  o{l—M2)c/  L  as  given  by  Poinsot  and 
Lele  in  [33]  for  Cartesian  coordinates.  The  corresponding  form  used  in  the  generalized 
coordinates  is 


K.  =  <j\\-  M2\/{V2JpL)  (5.27) 

where  M.  is  the  maximum  Mach  number  in  the  flow  field.  L  is  the  characteristic 
length  of  the  domain,  c  is  the  speed  of  sound.  The  preferred  range  for  constant  a  is 
0.2-0. 5.  The  absolute  value  of  1  —  M2  is  to  ensure  the  term  is  positive  because  the 
maximum  Mach  number  can  be  greater  than  1  in  a  transonic  flow  field. 

If  £5  =  0,  it  switches  to  the  “perfect”  non-reflective  boundary  condition.  However, 
this  boundary  condition  is  not  well  posed  and  will  not  lead  the  solution  to  the  one 
matching  the  exit  pressure  p^.  Equation  (5.26)  assumes  that  the  constant  exit  pres¬ 
sure  Poo  is  imposed  at  infinity.  There  exists  reflection  if  p  ^  p^,  which  is  needed  for 
the  well  posedness  of  the  numerical  solution.  For  the  unsteady  problems,  Equation 
(5.26)  will  make  the  mean  value  of  the  pressure  at  the  exit  very  close  to  p^.  However, 
the  pressure  at  the  individual  control  volume  may  not  be  exactly  equal  to  p^  even 
though  the  value  of  £5  can  be  very  small.  In  this  sense,  Equation  (5.26)  may  be 
considered  as  “almost  non-reflective  boundary  conditions” . 

The  complete  boundary  conditions  used  at  the  exit  are  the  pressure  at  infinity  for 
Equation  (5.26)  and  three  zero  gradient  viscous  conditions: 

(CzRcy  +  £>yTyy  +  £zrzy)  ~  0  (5.28) 

{ixTxz  +  ZyTyx  +  izTzz)  =  0  (5.29) 

^(txQx-HyQy  +  &QZ)=0  (5.30) 

The  amplitudes  of  the  outgoing  characteristic  waves,  £1,  £2,  £3,  and  £4  are  com¬ 
puted  from  the  interior  domain.  All  the  conservative  variables  at  phantom  points  are 
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obtained  by  solving  the  characteristic  N-S  equations,  Equation  (5.22).  All  the  trans¬ 
verse  and  viscous  terms  in  Equation  (5.22)  can  be  evaluated  in  the  same  way  as  the 
inner  domain  control  volumes.  The  Roe’s  Riemann  solver  is  also  used  for  computing 
fluxes  F'  and  G',  central  differencing  is  used  for  fluxes  E(,,  F(,,  F(,.  This  strategy 
makes  maximum  use  of  the  existing  code  and  minimizes  the  programming  work  in 
implementing  the  boundary  conditions. 

5.2.3  Subsonic  inflow  boundary  conditions 

At  £  =  1  boundary,  four  characteristic  waves,  £lt  £-2,  £3,  and  £4  are  entering  the 
domain  while  £5  is  leaving  the  domain.  For  3-D  open  field  flow  cases,  four  physi¬ 
cal  boundary  conditions  are  needed,  i.e.  u0 ,  v0,  wQ  and  p0  are  set  to  be  constant. 
Other  primitive  variables  are  specified  according  to  the  freestream  condition.  The 
total  energy  pe0  is  obtained  by  solving  the  energy  equation  in  Equation  (5.22).  The 
outgoing  wave  £5  can  be  estimated  by  using  interior  variables.  The  rest  of  the  waves 
are  evaluated  by  using  the  LODI  relations,  Equation  (5.23).  £\  -  £4  can  be  expressed 
as 


£\  —  —  lx-^~c{£i+ £^)i 


£2  — 


-?JL 

^yV2< 


(£4 +£5),  £3  —  -iz—^-(£i+ £5), 
V2c 


£4  —  £5 
(5.31) 


5.2.4  Adiabatic  wall  boundary  conditions 

At  a  3-D  adiabatic  wall  (77  =  constant),  the  no-slip  condition  is  enforced  by  extrapo¬ 
lating  the  velocity  between  the  phantom  and  interior  cells,  ua  =  —ut,  vQ  =  —  u,,  and 
wa  =  —Wi.  One  more  physical  boundary  condition  to  be  imposed  on  the  wall  is  the 
adiabatic  condition,  =  0.  From  the  adiabatic  condition,  the  p0  can  be  expressed 
as  the  following 


Po=Pi 

Po  Pi 


(5.32) 


The  total  energy  peQ  is  determined  by  solving  the  energy  equation  in  Equation 
(5.22).  Then  using  Equation  (5.32)  and  Equation  (??),  p0  and  p0  can  be  solved. 
Cross  a  rj  boundary,  vector  £  is  expressed  as  the  following: 
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v[$.£(|) + -  *£(J)  -  H(5»  <633) 

O'  +  OIMB)  +  ^^(5)  +  M<J)  +^(?)1 
(v  -  0[-$&(5)  -  $£(5)  -  $£(5) +/»£(!)!/ 

where  V'  =  ‘r].,u  +  qyv  +  //.Jr  and  C  =  cjVr)|,  V//  —  ,/rfe2  +  I],/  +  rj;2.  It  can  be  seen 
from  Equation  (5.33),  the  characteristic  waves  C\  -  £3  vanish  since  V  =  0  at  wall 
surface.  At  lower  wall  (t?  =  1),  the  outgoing  characteristic  wave  £5  is  computed  from 
the  interior  domain.  The  incoming  wave  £4  is  estimated  by  using  LODI  relations. 
By  solving  2nd  -  4th  equations  in  Equation  (5.23),  it  yields  £4  =  £5.  At  upper  wall 
(maximum  77),  the  £4  becomes  the  outgoing  wave,  and  it  can  be  computed  from  the 
interior  domain.  £5  is  the  incoming  wave  which  is  evaluated  by  £5  =  £4. 
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Chapter  6 

Structural  Models 

6.1  Modal  Approach  for  3D  Wing[4] 

The  governing  equation  of  the  solid  structure  motion  can  be  written  as, 

„  .cPu  du 

Mrfp'+C*'+Ku'f  <61) 

where  M,  C  and  K  are  the  mass,  damping,  and  stiffness  matrices  of  the  solid  respec¬ 
tively,  u  is  the  displacement  vector  and  f  is  the  force  exerted  on  the  surface  node 
points  of  the  solid,  both  can  be  expressed  as: 


/  Uj  > 

( fl  \ 

u  = 

Ui 

,f  = 

f i 

\  U  N  / 

V  f N  ) 

where  N  is  the  total  number  of  node  points  of  the  structural  model,  u*  and  f*  are 
vectors  with  3  components  in  x,  y,  z  directions: 


1  Uix  ^ 

l  1 

(f»\ 

Ui  = 

uiy 

II 

V  u-  ) 

'  1 

fi  is  dynamic  force  exerted  on  the  surface  of  the  solid  body.  In  a  modal  approach, 
the  modal  decomposition  of  the  structure  motion  can  be  expressed  as  follows: 
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K$  =  M$A 


(6.2) 


or 


K  <S>j  =  (6.3) 

where  A  is  eigenvalue  matrix,  A  =  diag[ Ai,  •  •  • ,  A A3at],  and  jth  eigenvalue  A j  = 
u>j2,  u>j  is  the  natural  frequency  of  jth  mode,  and  the  mode  shape  matrix  $  = 
[4> i> '  •  •  ‘ ' '  ,4>3 tv]. 

Equation  (6.15)  can  be  solved  by  using  a  finite  element  solver  (e.g.  ANSYS)  to 
obtain  its  finite  number  of  mode  shapes  <f)j.  The  first  five  mode  shapes  will  be  used 
in  this  paper  to  calculate  the  displacement  of  the  structure  such  that, 

“(<)  =  (6-4) 

3 

where  a  =  [a1,a2,a3,a4,a5]T.  Substitute  Equation  (6.4)  to  Equation  (6.1)  and  yield 

A  +  K$a  =  f  (6.5) 

dtz  dt 

Multiply  Equation  (6.5)  by  3>T  and  re-write  it  as 

*-clrdL  *  dci  . 

Mdp-+c*+Ka=p  (6-6) 

where  P  =  [Pi,  P2,  ■  •  • ,  Pj,  •  ■  • ,  Pn]t ,  the  modal  force  of  jth  mode,  Pj  —  <f>J f,  the 
modal  mass  matrix  is  defined  as 


M  =  =  diag{m\ ,  •  •  • ,  rrij,  •  •  • ,  m (6.7) 

where  rrij  is  the  modal  mass  of  jth  mode,  and  the  modal  damping  matrix  is  defined 
as 

C  =  $TC<f>  =  diag(ci ,  •  •  • ,  Cj,  ■  •  • ,  c3w)  (6.8) 

where  Cj  is  the  modal  damping  of  jth  mode,  and  the  modal  stiffness  matrix  is  defined 
as 


K  =  =  diag(ki,  k3N) 


(6.9) 
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where  kj  is  the  modal  stiffness  of  jth  mode.  Equation  (6.6)  implies 


(Paj 

~dt? 


dO-i  9 

+  2G  +UJjai 


(6.10) 


where  Q  is  modal  damping  ratio.  Equation  (6.18)  is  the  modal  equation  of  struc¬ 
ture  motion,  and  is  solved  numerically  within  each  iteration.  By  carefully  choosing 
reference  quantities,  the  normalized  equation  may  be  expressed  as 


where  the  dimensionless  quantities  are  denoted  by  an  asterisk,  uQ  is  the  natural 
frequency  in  pitch,  bs  is  the  streamwise  semichord  measured  at  wing  root,  L  is  the 
reference  length,  fh  is  the  measured  wing  panel  mass,  v*  is  the  volume  of  a  conical 
frustum  having  streamwise  root  chord  as  lower  base  diameter,  streamwise  tip  chord  as 
upper  base  diameter,  and  panel  span  as  height,  V*  =  and  Uoo  is  the  freestream 

velocity. 

Then  the  equations  are  transformed  to  a  state  form  and  expressed  as: 


[M]^>  +  [K]{S}  =  q 

where 

8-(2).M-[Z].K-((iJ).  4))-  = 


(6.12) 


To  couple  the  structural  equations  with  the  equations  of  flow  motion  and  solve  them 
implicitly  in  each  physical  time  step,  above  equations  are  discretized  and  integrated 
in  a  manner  consistent  with  Equation  (4.14)  to  yield 


(-U  +  +  k)  5Sn+l'm+1  =  _ M3Sn+1’m  ^Sn  +  S!..!  _  KSn+i,m  + 

VAr  At  )  2A  t 


n+l, m+l 


(6.13) 


where  n  is  the  physical  time  level  index  while  m  stands  for  the  pseudo  time  index. 
The  detailed  coupling  procedure  between  the  fluid  and  structural  systems  is  given  in 
the  following  chapter. 
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6.2  Mistuned  Bladed  Structural  Model  for  Tran¬ 
sient  Response 

The  Subsets  of  Nominal  Modes  (SNM)  structural  model  suggested  by  Yang  and  Grif¬ 
fin  [43]  is  developed  in  this  research  for  time  domain  to  calculate  the  structural  modes, 
which  are  expensive  to  calculate  if  the  direct  finite  element  approach  is  used.  Yang 
and  Griffin  recognized  that  each  mistuned  structural  mode  can  be  well  represented 
by  a  subset  of  the  tuned  structural  modes.  The  SNM  approach  was  developed  to 
take  the  finite  element  modal  solution  of  the  tuned  structure  as  the  input  to  formu¬ 
late  a  reduced  order  model  for  the  mistuned  structure.  The  order  of  the  problem 
thus  dropped  from  millions  to  hundreds,  and  the  computational  time  to  compute  a 
mistuned  structural  mode  is  reduced  from  hours  to  seconds.  This  is  critical  to  the 
simulation  of  fully  coupled  fluid-structural  problems  because  only  very  limited  com¬ 
putational  resource  is  required  in  addition  to  the  CPU  intensive  CFD  computation. 
A  brief  description  of  the  SNM  model  is  in  the  following: 

a)  Transformation  from  Finite  Element  Domain  to  Modal  Domain 

First,  the  equation  of  motion  in  the  finite  element  form  is  transformed  into  the 
modal  coordinates.  Assuming  that  the  variation  of  the  mechanical  damping  is  negli¬ 
gible,  then 


(M°  +  AM)5  +  C°a  +  (A0  +  AK)a  =  p  (6. 14) 

In  eq.  (6.14),  the  modal  coordinate  vector  a  is  the  displacements  of  the  tuned 
modes,  M°,  C°,  and  K°  are  the  modal  mass,  damping,  and  stiffness  matrices  of  the 
tuned  system,  and  typically  diagonal,  A K  and  AM  are  the  changes  in  modal  stiffness 
and  mass  matrices,  and  p  is  the  modal  force  vector.  Eq.  (6.14)  can  then  be  cast  in  a 
state  space  form, 


By  =  Ay  +  q  (6.15) 

where 

B  =  (  0  (M°  +  AM)  )  ’  y  =  (  a  )  (6'16) 

A  =  (  -(K°  +  AK)  -C0)’  q=(p)  ^6'17^ 

Eq.  (6.15)  is  the  modal  equation  of  motion  for  the  mistuned  structure.  Without 
truncating  the  modes,  the  order  of  eq.  (6.15)  is  2N  where  N  is  the  number  of  degrees  of 
freedom  of  the  whole  wheel  finite  element  model.  However,  the  order  can  be  reduced 
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to  2n,  where  n  is  tuned  modes  selected  in  the  SNM  representation,  to  simulate  the 
mistuned  structural  vibration  with  sufficient  accuracy.  N  (millions)  is  typically  much 
greater  than  n  (hundreds). 

b)  Diagonalization  of  the  Modal  Governing  Equation 

To  improve  the  computational  efficiency,  the  solution  of  eq.  (6.15)  can  be  further 
expressed  in  terms  of  its  right  eigenvectors  R  which  satisfies  the  following  eigenvalue 
problem 


AR  =  BRA  (6.18) 

where  A  is  the  diagonal  eigenvalue  matrix  of  the  mistuned  structure.  Applying  the 
classical  modal  analysis  technique,  eq.  (6.15)  can  be  transformed  in  a  diagonal  form 

B0  =  A0  +  q  (6.19) 

where  the  diagonal  matrices  B  and  A  are  the  generalized  mass  and  stiffness  ma¬ 
trices,  0  is  the  generalized  coordinates,  and  q  is  the  generalized  forces.  Since  the 
components  of  0  are  decoupled  from  each  other,  eq.  (6.19)  can  be  simulated  at  very 
low  computational  costs. 

Note  that,  in  eq.  (6.19),  B  and  A  mathematically  define  a  mistuned  bladed  disk 
structure,  the  generalized  force  q  is  derived  from  the  pressure  distribution  on  the 
airfoil  surfaces,  and  the  time-varying  unknown  0  will  be  solved  at  each  time  step. 
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Chapter  7 

Fully  Coupled  Fluid- Structural 
Interaction 


To  rigorously  simulate  fluid-structural  interactions,  the  equations  of  flow  motion  and 
structural  response  need  to  be  solved  simultaneously  within  each  iteration  in  a  fully 
coupled  numerical  model.  The  calculation  based  on  fully  coupled  iteration  is  CPU 
expensive,  especially  for  three  dimensional  applications.  The  modal  approach  can 
save  CPU  time  significantly  by  solving  the  modal  displacement  equations,  Eq.  (6.18), 
instead  of  the  original  structural  equations,  Eq.  (6.15),  which  is  usually  solved  by 
using  finite  element  method.  In  the  modal  approach,  the  structural  mode  shapes 
can  be  pre-determined  by  using  a  separate  finite  element  structural  solver.  Once 
the  several  mode  shapes  of  interest  are  obtained,  the  physical  displacements  can  be 
calculated  just  by  solving  those  simplified  linear  equations,  i.e.,  Eqs.  (6.18)  and 
(6.4).  In  present  study,  the  first  five  mode  shapes  provided  in  Ref. [55]  are  used 
to  model  the  wing  structure.  These  pre-calculated  mode  shapes  are  obtained  on  a 
fixed  structural  grid  system  and  are  transformed  to  the  CFD  grid  system  by  using  a 
3rd  order  polynomial  fitting  procedure.  The  procedure  is  only  performed  once  and 
then  the  mode  shapes  for  CFD  grid  system  are  stored  in  the  code  throughout  the 
simulation. 

The  procedure  of  the  fully  coupled  fluid-structure  interaction  by  modal  approach 
is  described  below: 

(1)  The  flow  solver  provides  dynamic  forces  on  solid  surfaces. 

(2)  Integrate  fluid  forces  at  each  surface  element  to  obtain  the  forcing  vector  f. 

(3)  Use  Eq.  (6.18)  to  calculate  modal  displacements  aj(j  =  1,2, 3, 4, 5)  of  the  next 
pseudo  time  step. 

(4)  Use  Eq.  (6.4)  to  calculate  physical  displacement  u  of  the  next  pseudo  time 
step. 
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(5)  Check  the  maximum  residuals  of  both  solutions  of  the  flow  and  the  structural 
equations.  If  the  maximum  residuals  are  greater  than  the  prescribed  convergence 
criteria,  go  back  to  step  (1)  and  proceed  to  the  next  pseudo  time.  Otherwise  the 
calculation  of  the  flow  field  and  the  structural  displacement  within  the  physical  time 
step  is  completed  and  the  next  new  physical  time  step  starts.  The  procedure  is  also 
illustrated  in  the  flow  chart  given  in  Fig.  12.132. 


Figure  7.1:  Fully  coupled  flow-structure  interaction  procedure 


Chapter  8 

Results  and  Discussion 


8.1  Validation  of  Zha-Hu  E-CUSP  Schemes [1] 

The  original  Zhu-Hu  E-CUSP  scheme  described  from  Eq.  4.25  to  4.35  was  developed 
and  then  used  as  the  basis  for  the  Zha  E-CUSP2  scheme. 

8.1.1  Shock  Tubes 

For  shock  tube  problems,  the  interests  are  focused  on:  1)  the  quality  (monotonicity 
and  sharpness)  of  the  shock  and  contact  discontinuities;  2)the  maximum  allowable 
CFL  number  to  be  used  for  explicit  Euler  method. 

For  explicit  Euler  time  marching  scheme,  it  is  desirable  that  the  CFL  number  is 
close  to  the  upper  limit  of  1.0.  For  the  ID  linear  wave  equation  with  CFL=1  and  1st 
order  upwind  scheme,  the  numerical  dissipation  and  dispersion  vanish.  For  nonlinear 
Euler  equations,  it  is  also  true  that  the  closer  the  CFL  to  1.0,  the  less  the  numerical 
dissipation. 


The  Sod  Problem 

Fig.  12.1  to  12.5  are  the  computed  temperature  distributions  using  different  upwind 
schemes  with  first  order  accuracy  compared  with  the  analytical  result  of  the  Sod 
problem  [56].  Since  the  computation  stops  before  the  waves  reach  either  end  of  the 
shock  tube,  the  first  order  extrapolation  boundary  conditions  are  used  at  both  ends 
of  the  shock  tube  for  all  the  schemes. 

The  maximum  allowable  CFL  number  for  a  scheme  is  defined  as:  beyond  which 
the  solution  will  either  be  oscillatory  or  unstable.  The  new  E-CUSP  scheme  (Zha 
CUSP  in  the  figures)  achieves  maximum  CFL  of  1.00,  and  the  shock  profile  is  the 
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crispest  and  remains  monotone  (Fig. 12.1).  The  maximum  allowable  CFL  of  Roe  and 
Van  Leer  scheme  are  0.95  and  0.96  respectively.  The  new  E-CUSP  scheme  takes  three 
grid  points  across  the  shock  wave,  while  the  Roe  and  Van  Leer  schemes  take  four  grid 
points  (see  Fig.12.1,  12.2  and  12.3).  The  Van  Leer  scheme  generates  a  tail  at  the  end 
of  the  expansion  wave  (see  Fig.  12.3).  Interestingly,  the  Van  Leer-Hanel  scheme  can 
reach  maximum  CFL  =1.0  and  the  shock  profile  is  also  crisper  than  the  original  Van 
Leer  scheme  with  no  tail  generated  at  the  end  of  the  expansion  wave  (see  Fig.  12.4). 
All  the  schemes  smear  the  contact  surface  to  a  similar  extent.  The  expansion  wave 
is  captured  well  by  all  the  schemes.  The  AUSM+  scheme  has  the  unexpectedly  low 
maximum  allowable  CFL  of  0.275.  The  whole  shock  and  contact  surface  profiles  are 
seriously  smeared  due  to  the  low  maximum  CFL  number. 

The  table  8.1  given  below  summarizes  the  maximum  allowable  CFL  number  for 
each  scheme.  Overall,  for  the  Sod  ID  shock  tube  problem,  the  new  scheme  suggested 
in  this  paper  performs  the  best  based  on  the  shock  sharpness,  monotonicity,  and 
stability. 


Table  8.1:  Maximum  CFL  Numbers  for  Sod  ID  Shock  Tube 
Scheme  CFL  Number 


The  new  scheme  (Zha  CUSP) 

1.00 

Van  Leer-Hanel 

1.00 

Van  Leer 

0.96 

Roe 

0.95 

Liou  AUSM+ 

0.275 

Slowly  Moving  Contact  Surface 

This  is  a  shock  tube  case  used  in  [57]  to  demonstrate  the  capability  of  the  scheme  to 
capture  the  contact  surface.  The  initial  conditions  are  [p,u,p]i  =  [0.125,0.112, 1.0], 
[p,  u,p\n  =  [10.0,0.112, 1.0]  .  All  the  results  are  first  order  accuracy.  Fig.  12.6  shows 
that  the  new  E-CUSP  scheme  ,  the  Roe  scheme  and  the  AUSM+  scheme  all  can  resolve 
the  contact  surface  accurately  as  they  are  designed.  The  results  of  those  schemes  are 
at  time  level  0.01.  The  velocity  is  uniformly  constant  and  the  density  discontinuity  is 
monotone.  The  new  E-CUSP  (Zha  CUSP)  scheme  has  far  higher  CFL  number  than 
the  other  schemes  with  the  value  of  1.00.  The  Roe  scheme  has  the  max  CFL=0.3, 
and  Liou’s  MUMS+  has  0.48.  Fig. 12.7  shows  that  the  Roe  scheme  generates  large 
velocity  oscillations  when  CFL=0.35,  greater  than  its  max  CFL=0.3. 

The  schemes  of  Van  Leer,  Van  Leer-Hanel  severely  distort  the  profiles  of  the  contact 
surfaces  as  shown  in  Fig.  12.8.  The  velocity  profiles  are  largely  oscillatory.  The 
density  jumps  are  also  more  smeared. 
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The  table  8.2  lists  the  maximum  CFL  number  of  each  scheme  for  the  slowing 
moving  contact  surface.  Again,  the  new  scheme  outperforms  the  other  schemes  by 
having  the  highest  CFL  number  and  still  maintain  the  monotonicity. 

Table  8.2:  Maximum  CFL  numbers  of  the  schemes  resolving  the  contact  surface 

Scheme  CFL  Number 

The  new  E-CUSP  (Zha  CUSP)  scheme  1.00 
Liou  AUSM+  0.48 

Roe  0.32 

Van  Leer  fail 

Van  Leer-Hanel  fail 


8.1.2  Entropy  condition 

This  case  is  to  test  if  a  scheme  violates  the  entropy  condition  by  allowing  the  expansion 
shocks.  The  test  case  is  a  simple  quasi-ID  converging-diverging  transonic  nozzle[58, 
59].  The  correct  solution  should  be  a  smooth  flow  from  subsonic  to  supersonic  with  no 
shock.  However,  for  an  upwind  scheme  which  does  not  satisfy  the  entropy  condition, 
an  expansion  shock  may  be  produced. 

For  the  subsonic  boundary  conditions  at  the  entrance,  the  velocity  is  extrapolated 
from  the  inner  domain  and  the  other  variables  are  determined  by  the  total  temper¬ 
ature  and  total  pressure.  For  supersonic  exit  boundary  conditions,  all  the  variables 
are  extrapolated  from  inside  of  the  nozzle.  The  analytical  solution  was  used  as  the 
initial  flow  field.  Explicit  Euler  time  marching  scheme  was  used  to  seek  the  steady 
state  solutions.  All  the  schemes  use  first  order  differencing. 

Fig.  12.9  is  the  comparison  of  the  analytical  and  computed  Mach  number  distribu¬ 
tions  with  201  mesh  points  using  the  new  scheme  and  the  scheme  of  Roe,  Van  Leer, 
Van  Leer-Hanel,  Liou’s  AUSM+.  The  analytical  solution  is  smooth  throughout  the 
nozzle  and  reaches  the  sonic  speed  at  the  throat  (the  minimum  area  of  the  nozzle, 
located  at  X/h  —  4.22).  It  is  seen  that  both  the  Roe  scheme  and  Van  Leer  scheme 
generate  a  strong  expansion  shock  at  the  nozzle  throat.  Both  schemes  can  converge 
to  machine  zero  (12  order  of  magnitude)  with  CFL=0.95  even  with  the  expansion 
shock  waves. 

The  Van  Leer-Hanel  scheme  can  not  converge  even  with  CFL=0.01.  The  result 
plotted  in  Fig.  12.9  is  the  one  before  it  diverges.  It  shows  an  expansion  shock  with 
the  Mach  number  jumping  from  0.74  to  1.42.  The  AUSM+  also  has  difficulties  to 
converge  for  this  case.  Using  CFL=0.05,  it  managed  to  reduce  the  residual  by  4  order 
of  magnitude.  The  solution  of  the  AUSM+  also  shows  an  expansion  shock  with  the 
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Mach  number  jumping  from  0.86  to  1.17. 

The  new  E-CUSP  scheme  does  not  have  an  expansion  shock  wave  at  the  sonic  point, 
but  is  not  smooth  due  to  the  discontinuity  of  the  first  derivative  of  the  pressure  at  the 
sonic  point.  This  is  shown  as  a  small  glitch  at  the  sonic  point  in  fig.  12.9.  The  glitch 
does  not  affect  the  scheme  to  converge  the  solution  to  machine  zero  with  CFL=0.95. 

As  indicated  in  [58,  59],  the  amplitude  of  the  expansion  shock  decreases  when  the 
mesh  is  refined.  When  the  2nd  order  schemes  with  the  MUSCL  differencing  are  used, 
all  the  expansion  shock  waves  as  well  as  the  glitch  of  the  new  scheme  at  the  sonic 
point  disappear.  Since  this  paper  is  to  compare  the  original  Riemann  solver  schemes, 
no  entropy  fix[60]  that  can  remove  the  expansion  shock  of  Roe  schemes  was  used. 

8.1.3  Wall  Boundary  Layer 

To  examine  the  numerical  dissipation  of  the  new  scheme,  a  laminar  supersonic  bound¬ 
ary  layer  on  an  adiabatic  flat  plate  is  calculated  using  first  order  accuracy.  The  in¬ 
coming  Mach  number  is  2.0.  The  Reynolds  number  based  on  the  length  of  the  flat 
plate  is  40000.  The  Prandtl  number  of  1.0  is  used  in  order  to  compare  the  numer¬ 
ical  solutions  with  the  analytical  solution.  The  baseline  mesh  size  is  81x61  in  the 
direction  along  the  plate  and  normal  to  the  plate  respectively. 

Fig.12.10  is  the  comparison  between  the  computed  velocity  profiles  and  the  Blasius 
solution.  The  solutions  of  the  new  scheme  (Zha  CUSP),  Roe  scheme,  and  AUSM+ 
scheme  agree  very  well  with  the  analytical  solution.  The  Van  Leer  scheme  significantly 
thickens  the  boundary  layer.  The  Van  Leer-  Hanel  scheme  does  not  improve  the 
velocity  profile. 

Fig.  12. 11  is  the  comparison  between  the  computed  temperature  profiles  and  the 
Blasius  solution.  Again,  the  new  scheme  (Zha  CUSP),  Roe  scheme,  and  AUSM+ 
scheme  accurately  predict  the  temperature  profiles  and  the  computed  solutions  ba¬ 
sically  go  through  the  analytical  solution.  Both  the  Van  Leer  scheme  and  the  Van 
Leer-  Hanel  scheme  significantly  thicken  the  thermal  boundary  layer  similarly  to  the 
velocity  profiles. 

Table  8.3  shows  the  wall  temperature  predicted  by  all  the  schemes  using  the  base¬ 
line  mesh  and  refined  mesh.  The  predicted  temperature  value  by  the  Van  Leer  scheme 
has  a  large  error.  The  Van  Leer-  Hanel  scheme  does  predict  the  wall  temperature 
accurately  even  though  the  overall  profile  is  nearly  as  poor  as  that  predicted  by  the 
Van  Leer  scheme.  The  new  scheme,  Roe  scheme  and  AUSM+  scheme  all  predict  the 
temperature  accurately. 

All  the  results  mentioned  above  are  converged  based  on  mesh  size.  The  wall 
temperatures  using  the  refined  mesh  of  161x121  are  also  given  in  table  8.3.  There  is 
little  difference  between  the  results  of  the  baseline  mesh  and  the  refined  mesh.  The 
refined  mesh  does  not  help  to  reduce  the  large  numerical  dissipation  of  the  Van  Leer 
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scheme.  When  the  2nd  order  schemes  are  used,  both  the  velocity  and  temperature 
profiles  of  the  Van  Leer  scheme  and  Van  Leer-  Hanel  are  improved  (not  shown). 


Scheme 

Twan ,  Mesh  80  x  60 

Twaii,  Mesh  160  x  120 

Blasius  analytical  solution 

1.8000 

1.8000 

The  new  E-CUSP  (Zha  CUSP)  scheme 

1.8025 

1.8018 

Roe  scheme 

1.8002 

1.7996 

Liou  AUSM+ 

1.8000 

1.8000 

Van  Leer 

1.8328 

1.8333 

Van  Leer-Hanel 

1.7970 

1.7996 

Table  8.3:  Computed  non-dimensional  wall  temperature  using  first  order  schemes 
with  the  baseline  mesh  and  refined  mesh 


8.1.4  Transonic  Converging-Diverging  Nozzle 

To  examine  the  performance  of  the  new  scheme  in  two-dimensional  flow  and  the 
capability  to  capture  the  shock  waves  which  do  not  align  with  the  mesh  lines,  a 
transonic  converging-diverging  nozzle  is  calculated  as  inviscid  flow.  The  nozzle  was 
designed  and  tested  at  NASA  and  was  named  as  Nozzle  Al[61].  Third  order  accuracy 
of  MUSCL  type  differencing  is  used  to  evaluate  the  inviscid  flux  with  the  Minmod 
limiter.  Fig.  12. 12  is  the  computed  Mach  number  contour  using  the  new  E-CUSP 
scheme  with  the  mesh  size  of  175  x  80.  In  the  axial  direction,  there  are  140  mesh 
points  distributed  downstream  of  the  nozzle  throat,  where  the  oblique  shock  waves 
are  located.  The  grid  is  clustered  near  the  wall.  For  clarity,  the  coarsened  mesh  is 
drawn  as  the  background  with  the  Mach  contours  to  show  the  relative  orientation  of 
the  shock  waves  and  the  mesh  lines.  The  nozzle  is  symmetric  about  the  centerline. 
Hence  only  upper  half  of  the  nozzle  is  calculated.  The  upper  boundary  uses  the 
slip  wall  boundary  conditions  and  the  lower  boundary  of  the  center  line  uses  the 
symmetric  boundary  conditions. 

As  indicated  by  the  wall  surface  isentropic  Mach  number  distribution  shown  in 
fig.  12. 13,  the  flow  is  subsonic  at  the  inlet  with  the  Mach  number  about  0.22  and  is 
accelerated  to  sonic  at  the  throat,  and  then  reaches  supersonic  with  Mach  number 
about  1.35  at  the  exit.  Fig.12.12  shows  that  right  after  throat,  an  expansion  fan 
emanates  from  the  wall  and  accelerates  the  flow  to  reach  the  peak  Mach  number 
about  1.5.  Due  to  the  sharp  throat  turning,  an  oblique  shock  appears  immediately 
downstream  of  the  expansion  fan  to  turn  the  flow  to  axial  direction.  The  two  oblique 
shocks  intersect  at  the  centerline,  go  through  each  other,  hit  the  wall  on  the  other 
side,  and  then  reflect  from  the  wall.  Such  shock  patten  is  repeated  to  the  exit  and 
the  shock  strength  is  weakened  with  the  flow  going  downstream.  Fig.  12.13  shows 
that  the  isentropic  Mach  number  distributions  predicted  by  the  new  CUSP  scheme 
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and  the  Roe  scheme  agree  fairly  well  with  the  experiment.  The  new  E-CUSP  scheme 
and  the  Roe  scheme  have  virtually  indistinguishable  results. 

The  mesh  refinement  study  indicates  that  the  mesh  resolution  in  the  axial  direction 
does  not  affect  the  shock  resolution  much.  The  axial  mesh  size  of  280  downstream  of 
the  throat  yields  only  slightly  better  shock  resolution  than  the  size  of  70.  However, 
the  mesh  size  in  the  vertical  direction  dramatically  changes  the  shock  resolution.  The 
mesh  size  of  80  in  the  vertical  direction  yields  much  better  resolution  than  the  mesh 
size  of  50.  This  can  be  seen  from  the  isentropic  Mach  number  in  fig.  12. 13,  which 
shows  that  the  mesh  size  of  175  x  80  generates  much  sharper  profiles  than  those  of 
the  mesh  175x50  for  the  first  and  second  shock  reflections. 

For  this  transonic  nozzle  with  the  mesh  size  175  x  80  on  an  Intel  Xeon  1.7Ghz 
processor,  the  CPU  time  per  time  step  per  node  to  calculate  the  inviscid  flux  is 
2.5871  x  10~6s  for  the  new  scheme,  which  is  about  25%  of  the  CPU  time  of  1.0284  x 
10~5s  used  for  the  Roe  scheme.  This  is  a  significant  CPU  time  reduction. 


8.1.5  Transonic  Inlet-Diffuser 

To  examine  the  performance  of  the  new  scheme  for  shock  wave/turbulent  boundary 
layer  interaction,  a  transonic  inlet-diffuser[62]  is  calculated  as  shown  by  the  Mach 
number  contours  in  fig.12.14,  which  has  the  exit  back  pressure  equal  to  0.83  times  of 
the  inlet  total  pressure.  The  Reynolds  number  based  on  the  throat  height  is  4.38  x  105. 
The  Baldwin-Lomax[48]  algebraic  turbulence  model  is  used.  Third  order  accuracy  of 
MUSCL-type  differencing  with  the  Minmod  limiter  is  used  for  the  inviscid  fluxes  and 
the  second  order  central  differencing  is  used  for  the  viscous  terms. 

A  normal  shock  is  located  downstream  of  the  throat  as  shown  in  fig.12.14.  No  flow 
separation  is  generated  under  this  back  pressure.  The  baseline  mesh  size  is  100  x  60. 
When  Ui  is  held  as  constant  and  the  mesh  is  refined  in  both  the  horizontal  and 
vertical  direction,  the  results  have  little  variation  and  are  converged  based  on  mesh 
size.  All  the  inlet-diffuser  results  presented  in  this  paper  are  from  the  mesh  size  of 
100  x  120.  The  mesh  in  the  horizontal  direction  is  clustered  around  the  shock  location 
to  better  resolve  the  shock  profile. 

Fig.  12.15  is  the  comparison  of  the  upper  wall  surface  pressure  between  the  experi¬ 
ment  and  the  computation.  The  agreement  is  very  good  except  that  the  computation 
predicts  the  shock  location  a  little  downstream  of  the  experimental  shock  location  and 
the  shock  strength  a  little  too  strong.  It  is  found  that  the  shock  profile  is  sensitive 
to  the  yjU  The  yf  value  of  2,2  x  10~4,7  x  10“6  are  tested.  The  smaller  yf  yields 
a  little  closer  shock  location  to  the  experiment.  The  results  shown  in  fig.12.14  and 
12.15  have  the  y+  value  of  2  x  10-4.  The  small  y*  effect  is  believed  due  to  the  first 
order  extrapolation  of  the  pressure  on  wall  surface  instead  of  the  requirement  of  the 
turbulence  modeling.  In  the  region  with  no  shock,  the  first  order  pressure  extrapola- 
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tion  on  the  wall  is  insensitive  to  the  distance  of  the  first  cell  to  the  wall,  while  in  the 
shock  region  it  is  sensitive  due  to  the  large  streamwise  gradient.  As  indicated  in  fig. 
12.15,  the  Roe  scheme  predicts  the  shock  location  slightly  closer  to  the  experiment 
than  the  new  CUSP  scheme. 

When  the  back  pressure  is  reduced  to  0.72  times  of  the  inlet  total  pressure.  The 
normal  shock  is  stronger  and  the  flow  separation  is  induced.  The  same  mesh  as  the 
previous  case  is  used  for  this  case.  Fig.  12. 17  is  the  predicted  pressure  distribution 
compared  with  the  experiment.  Both  the  new  CUSP  scheme  and  the  Roe  scheme 
predict  the  shock  location  accurately,  but  the  shock  strength  predicted  is  too  strong. 
However,  the  new  scheme  has  the  pressure  profile  in  the  separation  region  downstream 
of  the  shock  noticeably  closer  to  the  experiment  than  that  predicted  by  the  Roe 
scheme. 

It  should  be  pointed  out  that  the  turbulence  modeling  is  a  critical  factor  for  the 
prediction  accuracy  of  the  shock  wave/turbulent  boundary  layer  interaction.  Hence 
the  discrepancy  between  the  calculation  and  experiment  shown  above  is  only  partially 
attributed  to  the  different  discretization  schemes. 

Fig.12.17  is  the  pressure  contours  computed  using  p^t/Pt  =  0.72  with  the  new 
scheme,  Roe  Scheme,  and  Liou’s  AUSM+  scheme.  A  curved  A  shock  is  formed  due 
to  the  shock  wave/turbulent  boundary  layer  interaction.  The  shape  of  the  Mach 
contours  of  the  new  scheme  (Zha  CUSP)  and  the  Roe  scheme  are  very  much  alike. 
The  contours  computed  by  the  AUSM+  scheme  has  significant  oscillations  near  the 
wall. 


8.2  Validation  of  the  Zha  E-CUSP2  Scheme[3] 

The  Zha  E-CUSP2  scheme  defined  from  Eq.  4.25  to  4.39  was  validated  for  the  shock 
tube  and  laminar  boundary  layer.  The  results  are  as  good  as  the  original  Zha-Hu 
scheme  or  better.  More  importantly,  the  Zha  E-CUSP2  scheme  has  cured  the  tem¬ 
perature  oscillation  of  the  original  Zha-Hu  scheme  as  shown  below. 


8.2.1  Transonic  Converging-Diverging  Nozzle 

To  examine  the  performance  of  the  new  scheme  in  two-dimensional  flow  and  the 
capability  to  capture  the  shock  waves  which  do  not  align  with  the  mesh  lines,  a 
transonic  converging-diverging  nozzle  is  calculated  as  inviscid  flow.  The  nozzle  was 
designed  and  tested  at  NASA  and  was  named  as  Nozzle  Al[61].  Third  order  accuracy 
of  MUSCL  type  differencing  is  used  to  evaluate  the  inviscid  flux  with  no  limiter. 

Fig.  12. 18  is  the  computed  Mach  number  contours  using  the  original  Zha  CUSP 
scheme  and  the  Zha  CUSP2  scheme  with  the  mesh  size  of  175  x  80.  The  nozzle  is 
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symmetric  about  the  centerline.  Hence  only  upper  half  of  the  nozzle  is  calculated. 
The  upper  boundary  uses  the  slip  wall  boundary  conditions  and  the  lower  boundary 
of  the  center  line  uses  the  symmetric  boundary  conditions.  The  Mach  contour  lines 
computed  by  the  two  schemes  look  very  much  the  same.  However,  if  the  temperature 
contours  near  the  wall  is  zoomed  in,  it  can  be  seen  that  the  temperature  contours 
computed  by  the  Zha  CUSP  scheme  has  large  oscillations  as  shown  in  Fig.  12.19,  (a). 
The  temperature  oscillation  exist  along  the  whole  upper  wall  and  lower  symmetric 
boundary.  The  temperature  oscillations  are  removed  by  the  Zha  CUSP2  scheme  as 
shown  in  Fig.  12.19,  (b).  All  the  other  flow  variables  are  smooth  for  the  new  scheme. 

The  reason  that  the  original  Zha  CUSP  scheme  has  no  temperature  oscillation 
shown  in  Fig.  2  for  the  flat  plate  ,  but  has  oscillation  for  the  inviscid  nozzle  may  be  the 
following:  1)  The  laminar  Navier-Stokes  equations  provide  some  physical  dissipation 
to  smoothen  the  flow;  2)  The  boundary  layer  gradient  may  generate  some  numerical 
dissipation  to  smoothen  the  flow;  3)  The  1st  order  scheme  (piecewise  constant)  is 
used  for  the  flat  plate.  The  1st  order  scheme  is  monotone  and  has  higher  numerical 
dissipation  than  the  3rd  order  scheme  used  for  the  nozzle. 


8.3  2D  Flow  Induced  Vibration 

8.3.1  Stationary  Cylinder 

The  flow  past  a  stationary  cylinder  is  used  as  an  unsteady  flow  validation  case.  The 
baseline  mesh  dimensions  are  120x80  in  circumferential  and  radial  directions.  The 
far  field  boundary  is  located  20  diameters  away  from  the  center  of  the  cylinder.  The 
Reynolds  number  based  on  the  free-stream  condition  and  cylinder  diameter  is,  Re 
=  500.  The  laminar  Navier-Stokes  equations  will  be  solved  due  to  the  low  Reynolds 
number. 

The  computed  drag  and  lift  coefficients  are  shown  in  Figure  12.23.  As  shown  in 
the  figure,  the  lift  oscillates  at  certain  frequency  in  terms  of  the  Strouhal  number. 
The  drag  coefficient  oscillates  with  twice  that  frequency.  The  mesh  refinement  study 
and  computed  Strouhal  number,  drag,  lift  and  moment  coefficients  are  listed  in  Table 

8.4 


Table  8.4  shows  that  the  solution  is  converged  based  on  mesh  size.  The  computed 
lift  frequency  by  Zha-Hu  CUSP  scheme  agrees  well  with  the  experimental  results  of 
Roshko[66]  and  Goldstein [64],  and  is  closer  to  the  experimental  results  than  the  one 
computed  by  by  Alonso  et  al.  [65] ,  which  uses  more  grid  points. 
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Table  8.4:  Results  of  Mesh  Refinement  Study  and  comparison  with  the  experiments 


Mesh  Dimension 

StCl 

Q 

cd 

120x80 

I7y:;l 

Egg! 

0.2197 

±1.1810 

1.4529±0.1305 

200x120 

1 

■9 

0.2246 

±1.2267 

1.4840±0.1450 

(Roshko  1954[63]) 
(Goldstein  1938[64]) 
384x96  (Alonso  1995(65]) 

0.46735 

0.2075 

0.2066 

0.23313 

1.14946(Qmox) 

1.31523(CdQW) 

8.3.2  Vortex-Induced  Oscillating  Cylinder 
Structure  Model  of  Elastic  Cylinder 

For  the  computations  of  the  vortex-induced  oscillating  cylinder,  which  is  elastically 
supported  as  shown  in  Figure  ??  so  that  it  oscillates  only  in  the  direction  aligned 
with  or  normal  to  the  incoming  flow,  the  structural  dynamic  equations  which  govern 
the  motion  of  the  cylinder  are: 


mx  ±  Cxx  ±  Kxx  —  D  (8.1) 

my  ±  Cyy  ±  Kyy  =  L  (8.2) 

These  equations  are  solved  implicitly  together  with  the  equations  of  flow  motion, 
Equation  (4.14),  in  a  fully  coupled  manner.  In  Equation  (8.1),  x,  x,  and  x  represent 
the  dimensionless  horizontal  acceleration,  velocity  and  displacement  of  the  moving 
object  respectively.  Similarly,  y,  y,  and  y  in  Equation  (8.2)  represent  their  corre¬ 
sponding  ones  in  vertical  direction,  m ,  L,  and  D  are  the  mass,  lift,  and  drag  per  unit 
span  respectively,  Cx  and  Cy  are  the  damping  coefficients  in  horizontal  and  vertical 
directions,  Kx  and  Ky  are  the  spring  constants  in  horizontal  and  vertical  directions. 
In  present  study,  this  ’self-excited  oscillators’  is  designed  to  have  the  same  response 
in  both  direction,  i.e.  Cx  =  Cy  and  Kx  =  Ky. 

If  the  normalization  procedure  is  applied  to  Equations  (8.1)  and  (8.2)  by  using  the 
same  reference  scales  of  those  used  for  the  equations  of  flow  motion,  the  following 
nondimensional  equations  are  obtained 


x  ±  2£ 

D*+(l 

2V<r-  2  r 

/i57T 

(8.3) 

»+2c(; 

D*+(i 

r  1  V  = - Cl 

MJ 

(8.4) 

where  £  is  the  nondimensional  structural  damping  coefficient  calculated  by  £  = 
u  is  the  reduced  velocity  defined  by  u  =  b  is  radius  of  the  cylinder, 


Cx.V 


lyj  mKX'V 
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lj  =  the  mass  ratio,  fis  =  — ?lfe2,  Cd  and  Ci  are  the  drag  and  lift  force 

coefficients  respectively.  Then  the  equations  are  transformed  to  a  matrix  form  and 
expressed  by 


[M] 


d{s} 

dt 


+  [K]{S}  =  q 


(8.5) 


where 


S  = 


/  x 
x 

y 
\  v 


\ 

,  M  =  [I],  K  = 


/  0 

(1): 

0 

0 


o 

o 

(  0  \ 

0  0 

—cd 

0  -1 

,  q  = 

flsTT  a 

0 

I)'  2<(i)J 

To  couple  the  structural  equations  with  the  equations  of  flow  motion  and  solve  them 
implicitly  in  each  physical  time  step,  above  equations  are  discretized  and  integrated 
in  a  manner  consistent  with  Equation  (4.14)  to  yield 


a 


.  I+^M  +  K 
At  At 


)  SSn+1 


>m+1  =  _]V1 


3Sn+l,m  _  4gn  +  gn-1 

2A t 


_  j£gn+l,m  _|_  qn+  l,m+l 

(8.6) 


where  n  is  the  physical  time  level  index  while  m  stands  for  the  pseudo  time  index. 
The  detailed  coupling  procedure  between  the  fluid  and  structural  systems  is  given  in 
section  4. 


After  validating  the  stationary  cylinder  vortex  shedding  flow,  the  cylinder  is  re¬ 
leased  to  be  controlled  by  the  structure  model  as  shown  in  Figure  ??.  The  corre¬ 
sponding  structural  equations  are  given  in  section  3.1.  The  laminar  Navier-Stokes 
equations  are  solved  due  to  the  low  Reynolds  number,  Re  =  500. 

Using  the  temporally  periodic  solution  obtained  in  the  computation  of  stationary 
cylinder  as  the  initial  flow  field,  the  computation  is  resumed  after  the  cylinder  is  let 
to  move  in  both  streamwise  and  transverse  directions.  For  the  purpose  of  comparison 
with  the  experimental  data  of  Griffin  [67]  several  different  combinations  of  structural 
parameters  are  used  in  the  computations. 

Morton  et  al.[30]  suggested  to  use  the  reduced  velocity  u  —  such  that  the 
structural  oscillator  works  under  or  near  the  resonance  conditions.  Therefore  the 
computed  St  number  from  the  stationary  cylinder  is  used  to  determine  u.  For  all  the 
cases  of  oscillating  cylinder,  St  is  set  to  be  0.2,  corresponding  to  u  =  1.5915.  Different 
mass  ratios,  /is,  are  used  to  test  the  different  responses  of  the  structural  system.  They 
are  equal  to  1.2732,  5.0,  and  12.7324  respectively.  To  match  the  wide  range  of  the 
experimental  data,  the  damping  ratio,  (,  is  chosen  from  the  range  between  0.001  - 
1.583. 


8.3.  2D  FLOW  INDUCED  VIBRATION 


59 


The  dimensionless  physical  time  step  A t  —  0.05  is  used,  which  corresponds  to 
approximately  100  time  steps  per  period  determined  by  the  Strouhal  number  used. 
The  CFL  number  for  the  pseudo  time  steps  varies  from  case  to  case.  For  the  large 
cylinder  movement  cases,  smaller  pseudo  time  steps  are  used  to  limit  the  displacement 
of  the  cylinder  during  each  iteration. 

For  the  cases  computed,  the  CFL  number  varies  from  5  to  500.  The  iteration 
number  within  one  physical  time  step  varies  from  20  to  100.  Fig.  12.24  shows  a 
typical  iteration  history  within  one  physical  time  step.  Both  the  residual  of  the 
CFD  solver  and  structure  are  reduced  to  machine  zero.  The  structure  solver  usually 
converges  faster  than  the  CFD  solver. 

Figure  12.25  displays  the  computed  vorticity  contours  around  the  oscillating  cylin¬ 
der.  It  shows  how  the  vortexes  are  shed  at  the  moment  when  the  cylinder  bounds 
back  toward  its  mean  position  in  y  direction. 

A  typical  trajectory  of  the  center  position  of  the  moving  cylinder  is  plotted  in 
Figure  12.26,  which  is  similar  to  the  results  computed  by  Blackburn  et  al. [68]  and 
Alonso  et  al. [65] 

All  the  numerical  results  for  present  study  are  plotted  in  Figure  12.27  for  the  three 
values  of  [is.  Also  plotted  are  the  computations  of  Alonso  et  al.  [65]  with  ys  = 
5.0,  computations  of  Morton  et  al.  [30]  with  y,s  —  12.73,  and  the  experimental  data 
of  Griffin[67].  In  the  figure,  the  abscissa  is  the  reduced  damping  with  the  form  of 
87r25t2C~§2  [68],  and  the  ordinate  is  the  cross-flow  displacement  of  motion  normalized 
by  the  cfiameter  of  the  cylinder.  Overall,  very  good  agreement  is  observed  between 
the  present  results  and  the  experimental  results,  especially  for  the  case  of  y,s  =  1.2732. 
The  figure  shows  that  the  higher  values  of  mass  ratios  {jis  =  5.0  and  y,s  =  12.7324) 
give  less  satisfactory  results  than  those  with  ys  =  1.2732,  particularly  at  low  damping 
ratios.  However,  they  agree  well  with  the  results  of  Morton  et  al. [30]  (p,5  =  12.73). 

8.3.3  Elastically  Mounted  Airfoil 

As  the  validation  of  the  Zha-Hu  CUSP  scheme  for  transonic  airfoils,  the  steady  state 
solution  of  the  transonic  RAE  2822  airfoil  is  calculated  first.  The  freestream  condition 
for  this  study  are  listed  in  Table  8.5  below. 


Table  8.5:  Free-stream  condition  for  RAE  2822  Airfoil 


Mach  number 

Static  Pressure  (psia) 

Temperature  (R) 

Angle-of-Attack  (deg) 

Reynolds  Number 

0.729 

15.8073 

460.0 

2.31 
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The  turbulent  Reynolds  stress  and  heat  flux  is  calculated  by  the  Baldwin-Lomax 
algebraic  model  [48].  This  case  is  run  using  an  O-type  grid  with  three  different  di¬ 
mensions,  they  are  128x50x1,  256x55x1,  and  512x95x1  respectively.  The  outer 
boundary  extends  to  15  chords  from  the  center  of  the  airfoil.  The  experimental  data 
provided  by  Cook  et  al.  [69]  are  available  for  validation.  The  comparison  of  pressure 
coefficient  on  the  airfoil  is  shown  in  Figure  12.28.  Overall,  very  good  agreement  is 
obtained  between  the  computation  and  experiment  for  each  mesh  dimension,  espe¬ 
cially  for  the  two  larger  ones  which  appear  to  be  sufficient  to  capture  the  shock  on  the 
suction  surface  of  the  airfoil  without  using  any  limiter.  The  important  aerodynamic 
coefficients  from  both  simulation  and  experiment  are  summarized  in  Table  8.6. 


Table  8.6:  Aerodynamic  coefficients  and  y+  for  RAE  2822  Airfoil 


Mesh  Dimension 

0* 

a 

cm 

y+ 

128x50 

0.01482 

0.73991 

0.09914 

0.0833  -  2.3864 

256x55 

0.01455 

0.73729 

0.09840 

0.1318  -  2.4016 

512x95 

0.01426 

0.74791 

0.09994 

0.2309  -  2.0228 

Prananta  et  al.[70] 

0.01500 

0.74800 

0.09800 

Experiment 

0.01270 

0.74300 

0.09500 

It  can  be  seen  in  Table  8.6  that  the  predicted  lift  coefficients  with  all  mesh  dimensions 
agree  well  with  the  experiment.  The  computed  drag  and  moment  coefficients  show 
larger  errors,  but  they  have  the  similar  accuracy  as  those  computed  by  Prananta  et 
al.[70].  The  relatively  large  error  of  the  drag  and  moment  may  be  mostly  due  to  the 
inadequacy  of  the  turbulence  model,  which  is  difficult  to  predict  the  surface  friction 
accurately. 


8.3.4  Forced  Pitching  Airfoil 

As  a  validation  case  of  the  scheme  for  moving  grid  system,  the  forced  pitching  NACA 
64A010  airfoil  is  calculated.  For  this  transonic  airfoil,  the  Reynolds  averaged  Navier- 
Stokes  equations  with  Baldwin-Lomax  turbulence  model  are  solved.  Similar  to  the 
previous  computation  of  the  flow  over  the  stationary  airfoil,  an  O-type  mesh  consisting 
of  280x65  cells  is  employed  for  the  computations  of  forced  pitching  airfoil.  The 
NACA  64A0101  airfoil  is  selected  for  this  calculation  because  the  experimental  data 
is  available.  The  fine  mesh  zone  or  the  non-deforming  part  of  the  mesh  is  shown  in 
Figure  12.29.  The  first  grid  point  adjacent  to  the  wall  has  the  maximum  y+  <  3.43 

The  NACA  64A0101  airfoil  is  forced  in  pitch  around  its  quarter  chord  sinusoidally. 
The  angle  of  attack  is  imposed  as  a  function  of  time  as  a(t)  =  am  +  a0sin(u>t),  where 
am  and  aQ  are  the  mean  angle  of  attack  and  the  amplitude  of  oscillation  respectively. 
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The  u>  is  the  angular  frequency  which  is  directly  related  to  the  reduced  frequency 
Kc  —  2^->  where  c  is  the  airfoil  chord,  and  is  the  free-stream  velocity.  To 
compare  with  the  experimental  results  given  by  Davis[71],  the  primary  parameters 
used  in  the  computation  are  listed  as  follows:  am  =  0,  aQ  =  1.01°,  Re  =  1.256  x  107, 
Moo  =  0.8,  reduced  frequency,  Kc  =  0.202. 

Again,  the  computation  begins  with  the  steady  state  flow  field  of  the  stationary 
airfoil  at  0  degree  angle  of  attack  with  a  dimensionless  time  step  At  =  0.1.  The 
transition  period  takes  about  3  cycles  and  the  results  becomes  periodic  in  time  after 
that.  Figure  12.30  shows  the  lift  oscillation  versus  the  angle  of  attack  after  the  flow 
field  reaches  its  temporally  periodic  solution.  The  computed  lift  oscillation  agrees 
well  with  the  experiment  [71],  which  has  an  evident  improvement  compared  to  the 
recent  result  computed  by  McMullen  et  al.  in  2002  [72]. 

Fig.  12.31  shows  the  computed  moment  coefficient  compared  with  the  experimental]. 
The  computed  moment  coefficient  does  not  agree  as  accurately  with  the  experiment 
as  the  lift  coefficient  does.  However  the  results  are  very  similar  to  those  predicted  by 
Bohbot  et  al.  [73]  and  McMullen  et  al.  [72].  The  large  discrepancy  between  the  compu¬ 
tation  and  the  experiment  for  the  moment  coefficient  may  be  due  to  the  inadequacy 
of  the  turbulence  modeling,  which  may  not  predict  the  surface  friction  accurately. 


8.3.5  Flow-Induced  Vibration  of  NACA  64A010  Airfoil 
Structural  Models 

The  structural  model  for  the  flow-induced  vibration  of  a  2-D  sweptback  wing  with  a 
NACA  64A010  cross-section  is  described  in  section.  This  model  was  first  introduced 
by  Isogai  [74]  [74],  and  has  been  numerically  investigated  by  several  researchers  [10] 
[75]  [70]  [73], 

The  system  of  the  elastically  mounted  airfoil  is  assumed  to  be  undamped.  The 
airfoil  is  allowed  to  move  in  pitch  about  a  given  elastic  axis  and  plunge  vertically.  The 
pitch  axis  is  defined  by  a  distance  a,  which  is  the  multiple  of  the  semi-chord  length 
with  the  origin  point  located  at  the  mid-chord  position.  If  a  is  positive,  it  means  the 
axis  is  located  downstream  of  the  mid-chord,  negative  means  being  located  upstream 
of  the  mid-chord  point. 

A  sketch  of  the  elastically  mounted  airfoil  is  depicted  in  Figure  12.21.  The  motion 
of  such  an  elastic  system  can  be  described  by  using  the  following  equations 

mh  +  Saa  +  Khh  =  —  L  (8.7) 

SQh  T  IQOi  T  AQa  =  M  (8-8) 

where  h  and  a  are  the  plunging  and  pitching  displacements  respectively,  m  is  the  mass 
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per  unit  span,  SQ  is  the  static  moment  around  the  elastic  axis,  Ia  is  the  rotational 
moment  of  inertia,  Kh  and  Ka  are  plunging  and  pitching  spring  constants  respectively, 
L  is  the  lift  force  and  M  is  the  moment  around  the  elastic  axis.  The  equations  of  the 
structure  motion  (8.7)  and  (8.8)  are  normalized  by  using  semi-chord  b  as  the  length 
dimension,  the  uncoupled  natural  frequency  in  pitch  ua  as  the  time  scale,  and  are 
expressed  as 


-•  .  (Uh\2  h  U*2r 

h  +  xaa+[ —  I  h  = - Ci 

\u >aJ  fiir 

U*2 

xji  +  rla  +  rla  = - Cm 

fl7f 


(8.9) 

(8.10) 


where  xa  is  the  static  unbalance,  ujh  is  the  uncoupled  natural  frequency  in  plunge,  r2a 
is  the  squared  radius  of  gyration,  U*  is  the  reduced  velocity  defined  as  Ci  and  Crn 
are  the  lift  and  moment  coefficient  respectively.  Since  the  time  scale  used  in  Equations 
(8.9)  and  (8.10)  is  different  from  the  one  used  in  the  governing  equations  of  flow,  the 
structural  dimensionless  time  t*  needs  to  be  re-scaled  and  keep  its  consistency  with  the 
entire  system  during  the  computation,  i.e.,  t*s  =  where  t*j  is  the  dimensionless 

time  for  flow  and  the  L  is  the  length  scale.  Finally  the  equations  are  cast  into  the 
form  of  Equations  (8.5)  and  (8.6),  and  the  corresponding  matrices  are 
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The  structural  parameters  used  in  this  model  are  listed  as  the  following:  a  =  —2.0, 
xQ  =  1.8,  ^  =  1,  =  3.48,  and  \i  =  60.  The  elastic  axis  is  located  half  a  chord 

upstream  of  the  airfoil  nose. 


The  unsteady  Reynolds  averaged  Navier-Stokes  equations  with  the  Baldwin-Lomax 
turbulence  model  are  solved  for  the  flow  field  in  this  study.  The  freestream  conditions 
are:  Re  =  1.256  x  107,  =  0.825. 


Due  to  the  symmetric  profile  of  the  NACA  64A010  airfoil,  an  initial  perturbation  is 
imposed  to  trigger  the  oscillating  motion.  The  airfoil  is  forced  to  rotate  sinusoidally 
about  its  elastic  axis  at  the  natural  frequency  in  pitch  uQ  with  an  angle  of  attack 
amplitude,  a0  =  1°.  Usually  the  forced  pitching  mode  lasts  for  1  -  3  cycles.  After 
that,  the  elastically  mounted  airfoil  is  let  to  move  in  both  plunging  and  pitching 
directions,  and  then  the  dynamic  response  is  recorded. 


In  present  study,  the  search  of  the  critical  point  on  the  transonic  flutter  boundary 
at  a  given  Mach  number  is  conducted.  The  speed  index,  V*  defined  as  is  the 
parameter  to  classify  damped,  neutral  and  divergent  responses  of  the  airfoil  when  the 
Mach  is  fixed.  In  this  case,  the  total  pressure  and  temperature  need  to  be  adjusted 
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to  match  the  certain  value  of  the  Re  number.  Several  calculations  are  needed  to 
determine  the  critical  point  using  a  bi-section  method. 

In  Figures  (12.34)  through  (12.36)  the  time  histories  of  plunging  and  pitching 
displacements  at  —  0.825  are  plotted  for  three  different  V*.  In  these  figures, 
from  V*  =  0.55  to  V*  =  0.70,  the  plots  correspond  to  the  damped,  neutral,  and 
diverging  responses  respectively.  The  major  task  of  calculating  a  flutter  boundary 
is  to  locate  where  the  neutrally  stable  (critical  point)  is  by  looking  at  those  plots 
and  determining  where  the  neutral  response  occurs  as  the  V*  varies.  When  the 
value  of  V*  is  smaller  than  the  critical  value  on  the  flutter  boundary,  both  plunging 
and  pitching  displacements  decay  corresponding  to  the  damped  response  as  shown  in 
Figure  (12.34).  Once  the  value  of  V*  coincides  with  or  is  close  to  the  critical  value, 
the  neutral  response  appears  as  shown  in  Figure  (12.35).  Any  value  of  V*  beyond 
the  critical  value,  a  diverging  response  is  expected  as  shown  in  Figure  (12.36).  Mach 
number  0.825  is  located  at  the  bottom  of  the  sonic  dip  as  reported  in[5,  70,  73].  The 
predicted  critical  velocity  index  V*  =  0.615  is  consistent  with  the  results  computed 
by  those  researchers. 

8.3.6  2D  Airfoil  Flutter  Boundary  Prediction 

After  the  validation  of  the  flow  induced  vibration  of  NACA  64A010  Airfoil  in  chapter 
8.3.5,  a  search  of  the  flutter  boundary  is  conducted.  This  is  ultimately  the  most 
important  results  needed  for  an  aircraft  and  engine  design  engineer  to  determine  if 
the  design  is  located  within  the  safety  margin. 

At  each  Mach  number,  several  calculations  are  needed  to  determine  the  critical 
point  on  the  flutter  boundary  using  a  bi-section  method.  At  certain  Mach  number, 
the  flutter  boundary  is  very  ‘thin’,  and  more  calculations  are  necessary  to  really 
capture  the  critical  points.  The  dynamic  response  immediately  after  the  transonic 
dip  becomes  very  complex,  and  locating  the  flutter  boundary  in  that  region  ( Mach 
—  0.875  -  0.9)  is  very  difficult  and  time-consuming. 

The  V*  and  the  frequency  ratio  ~  for  the  flutter  boundary  are  plotted  versus 
Mach  number  in  Figures  (12.37)  and  fl2.38)  respectively.  Also  plotted  in  the  figures 
are  the  results  from  two  other  computations  by  Prananta  et  al.  [70]  and  Bohbot  et 
al.  [73].  The  Mach  number  for  the  bottom  of  the  transonic  dip  of  0.825  is  consistent 
with  their  results. 

Overall,  the  present  results  compare  well  with  the  results  of  of  Prananta  et  al.  and 
Bohbot  et  al.  except  that  both  values  of  V*  and  ^  ratio  are  higher  at  high  Mach 
number  region  ( Mach  =  0.925  to  0.95).  The  primary  difference  between  the  present 
results  and  their  results  are:  1)  The  present  results  are  based  on  fully  coupled  fluid- 
structure  interaction.  Their  results  are  loosely  coupled;  2)  The  Reynolds  number  of 
the  present  results  is  about  twice  higher.  Both  the  present  results  and  Prananta  et 
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al.  are  based  on  the  same  Baldwin-Lomax  turbulence  model.  To  study  the  effect  of 
Reynolds  number,  the  flutter  boundary  at  Mach  =  0.925  and  0.95  are  re-calculated 
with  the  same  Reynolds  number  (6  x  106)  used  by  Prananta  et  al.  The  difference  is 
small.  Hence,  the  difference  between  the  present  results  and  the  results  of  of  Prananta 
et  al.  and  Bohbot  et  al.  may  be  due  to  the  fully  coupled  and  loosely  coupled  algorithm. 
Since  there  are  no  experimental  results  for  comparison,  it  is  difficult  to  judge  which 
result  is  more  correct. 

Except  the  flutter  boundary,  the  present  solver  also  has  captured  the  Limit  Cycle 
Oscillation  (LCO)  phenomenon  as  shown  in  Figure  (12.39).  LCO  occurs  when  the 
velocity  index  is  greater  than  the  flutter  velocity  index[76].  The  amplitude  is  large, 
but  stable.  LCO  is  considered  due  to  the  nonlinear  nature  of  the  shock  boundary  layer 
interaction  occurring  for  transonic  airfoil[76,  73].  Figure  (12.40)  presents  a  different 
LCO  with  the  second  torsion  mode  more  dominant. 

Figure  (12.41)  shows  an  interesting  situation  happening  at  Mach  =  0.875  and 
in  the  area  of  V*  =  2.5.  Both  plunging  and  pitching  displacements  of  the  system 
increase  rapidly  immediately  after  the  airfoil  is  set  to  be  free,  but  then  gradually 
reach  their  steady  state  positions  and  stay  there  through  the  end  of  the  computation. 
Under  this  flow  condition,  the  aerodynamic  forces  and  moments  are  balanced  by  the 
structure  system.  The  angle  of  attack  is  stabilized  at  2.9°  and  is  in  the  range  of  the 
cruise  point  which  should  be  stable.  This  situation  is  only  observed  at  Mach  =  0.875, 
and  maybe  named  as  ‘standing’.  The  V*  of  the  standing  phenomenon  is  not  located 
at  one  point,  but  a  region  around  V*  =  2.5.  The  results  and  flow  conditions  of  the 
‘standing’  phenomenon  needs  to  be  confirmed  by  experiment. 


8.4  SNM  Model  Used  for  Transient  Response 

A  structural  model  is  developed  to  simulate  the  transient  response  of  a  blade  disk 
structure.  This  model  implements  the  original  SNM  theory  (for  steady  state  re¬ 
sponses)  and  expands  it  to  simulate  transient  responses.  The  study  discussed  in  this 
section  assumes  that  the  fluid  forcing  function  is  prescribed  such  that  the  transient 
response  can  be  carried  out  in  the  transient  simulation  of  ANSYS.  The  prescribed 
forcing  function  however  does  have  multiple  frequency  components  with  all  but  one 
primary  frequency  components  decaying  in  time. 

Fig.  12.42  shows  a  finite  element  model  of  24-blade  disk  with  dynamic  character¬ 
istics  similar  to  a  typical  high  compressor  rotor.  The  finite  element  mesh  is  rather 
coarse  for  the  purpose  of  tractable  computational  cost  in  ANSYS  benchmark  com¬ 
putation.  When  studying  a  realistic  case,  the  mesh  needs  to  be  finer.  The  transient 
responses  of  this  model  were  computed  by  two  methods:  In  the  first  method,  the  tran¬ 
sient  responses  were  calculated  using  the  ’’Full  Solution”  function  in  ANSYS  which 
utilizes  full  system  matrices  (instead  of  reduced  modal  matrices  in  ’’Reduced  Solu- 
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tion”  function)  to  avoid  modal  truncation  error.  The  second  method  computes  the 
transient  responses  based  on  the  SNM  method.  The  SNM  results  are  then  compared 
with  the  ANSYS  benchmarks. 

Both  low  frequency  response  and  high  frequency  response  of  a  tuned  bladed  disk 
were  studied. 

8.4.1  Low  Frequency  Case 

The  blade  tips  are  driven  by  prescribed  forces  with  a  circumferential  pattern  of  zero 
nodal  diameter.  The  force  is  given  by  Eq.8.11, 


f(t)  =  fdc  +  fosin(u>t)  +  ficos(l.4u)t)eoji^wo°  +  f2sin(0.7ut)eu/2OO°  (8.11) 

with  fdc  —  0,/o  =  1,  /i  =  0.5,  f2  =  2,  and  t o  —  500Hz.  Note  that  the  lowest  natural 
frequency  of  all  zero  nodal  diameter  modes  is  495  Hz.  The  overall  time  history  of 
the  force  is  shown  in  Fig.  12.43.  Fig.12.44  shows  that  multiple  frequency  contents 
(0.7a;,  a ;,  and  1.4a;)  co-exist  in  the  forces  at  the  beginning  of  the  simulation  and  Fig. 
12.45  shows  that  only  the  primary  frequency  component  (a;)  survives  at  the  end  the 
simulation.  The  forcing  function  is  chosen  so  to  mimic  the  dynamic  pressure  of  a 
fluid  field  and  to  examine  the  accuracy  of  the  SNM  model. 

Fig.  12.46  shows  the  envelopes  of  the  responses  of  a  tip  node  of  a  blade  computed 
by  SNM  and  ANSYS.  24  modes  of  one  blade  mode  family  was  used  in  the  SNM 
simulation.  Closer  inspections  of  the  SNM  and  ANSYS  results  shown  in  Fig.  12.47 
and  Fig.  12.48  verify  that  SNM  is  very  accurate  in  predicting  the  transient  response 
for  this  case. 

A  worthy  note  is  that  the  numbers  of  degrees  of  freedom  in  the  ANSYS  and  SNM 
simulations  are  2256  and  24  respectively.  The  time  saving  of  SNM  is  about  2  orders 
of  magnitude  over  this  simplified  ANSYS  model.  In  the  case  of  a  realistic  ANSYS 
model,  the  time  saving  would  be  much  greater. 


8.4.2  High  Frequency  Case 

In  the  high  frequency  case,  the  blade  tips  are  again  driven  by  prescribed  the  forces 
defined  by  Eq.  (8.11)  with  a  circumferential  pattern  of  zero  nodal  diameter.  The 
parameters  of  the  force  are  fdc  =  0 ,/0  =  l,/i  :  0.5,  f2  —  2,  and  a;  =  2500  Hz. 
Note  that  the  3rd  natural  frequency  of  all  zero  nodal  diameter  modes  is  2506  Hz. 
The  overall  time  history  of  the  force  is  shown  in  Fig.  12.49.  Fig.  12.50  and  12.51 
show  that  multiple  frequency  contents  (0.7a;,  u>,  and  1.4a;)  of  the  force  co-exist  at  the 
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beginning  of  the  simulation  and  only  the  primary  frequency  component  (lo)  survives 
at  the  end  the  simulation. 

Fig.  12.52  compares  the  envelopes  of  the  responses  of  a  tip  node  of  a  blade  com¬ 
puted  by  SNM  and  ANSYS.  Only  24  modes  of  the  3rd  blade  mode  family  were  used 
in  the  SNM  simulation.  Fig.  12.53  and  Fig.  12.54  are  closer  examination  of  the 
responses.  It  is  apparent  that  SNM  does  not  capture  the  response  well  initially  when 
the  forcing  frequency  components  and  are  yet  to  die  out  but  does  capture  the  re¬ 
sponse  quite  well  at  the  end  when  the  primary  frequency  component  dominates.  This 
is  because  of  the  lack  of  low  and  high  frequency  modes  in  the  SNM  model.  In  this 
case,  the  SNM  model  actually  acts  as  a  frequency  filter  of  the  overall  system.  When 
long  term  steady  state  solution  is  of  the  only  interest,  an  SNM  model  with  a  single 
blade  mode  family  might  be  desirable. 

To  further  verify  the  above  statement,  an  SNM  model  with  the  3rd  blade  mode 
family  and  the  neighboring  families  (2nd  and  4th  blade  modes)  was  used  for  the 
computation.  Fig.  12.55,  12.56,  and  12.57  compare  the  responses  computed  by  SNM 
and  ANSYS.  They  show  good  agreement  between  SNM  and  the  ANSYS  benchmark. 
Notice  that  the  number  of  degrees  of  freedom  in  the  3-family  SNM  model  is  72.  The 
time  saving  of  SNM  is  about  1.5  orders  of  magnitude  over  the  ANSYS  model. 


8.5  Non-Reflective  Boundary  Conditions 

8.5.1  A  Vortex  propagating  through  a  outflow  boundary 

The  first  test  case  is  a  subsonic  vortex  propagating  flow  in  an  open  flow  field.  The 
computed  domain  is  rectangular  inclined  30°  about  the  horizontal  axis  as  a  validation 
for  the  generalized  coordinates.  The  subsonic  inflow  and  outflow  boundary  conditions 
are  used  at  the  inlet  and  exit,  and  far  filed  boundaries  are  used  on  the  upper  and  lower 
borders.  In  [33],  a  supersonic  vortex  propagating  flow  is  chosen.  It  is  known  that  a 
subsonic  vortex  propagating  flow  is  more  difficult  to  deal  with  since  the  disturbance 
propagates  both  upstream  and  downstream.  To  test  present  method  under  more 
general  conditions,  the  subsonic  vortex  propagating  flow  is  selected  for  this  study. 

The  computational  mesh  has  a  length  of  2  units  in  streamwise  direction  which  is 
30  degrees  to  the  x  direction,  while  the  width  is  4  units  in  transverse  direction.  The 
mesh  dimensions  are  60x100.  The  laminar  Navier-Stokes  equations  is  solved  for  this 
flow  with  M  =  0.8  and  Reynolds  number  of  300. 

A  vorticity  is  initially  located  at  the  center  of  the  domain  when  dimensionless  time 
t*  =  0,  and  convected  downstream  toward  the  outflow  boundary.  The  velocity  flow 
field  is  initially  specified  as 
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u 

v 


x2  +  y2\ 


(8.12) 

(8.13) 


where  u ^  and  vx  are  the  velocity  components  of  the  incoming  flow,  Cv  is  the  coef¬ 
ficient  that  determines  the  vortex  strength  of  the  velocity  field,  and  Rc  is  the  vortex 
radius.  The  total  energy  field  is  initialized  as 


<’e  =  pe~+$'exp{-!L2jif)  <8'14) 

where  Ce  is  the  coefficient  that  determines  the  vortex  strength  of  the  total  energy  field. 
Equations  (8.12)  and  (8.14)  are  adopted  from  those  used  by  Poinsot  and  Lele[33].  The 
coefficients  Cv  and  Ce,  and  radius  Rc  are  defined  by 


Cv/(cL)  =  —0.0005,  Ce/(cL)  =  - 0.02,  Rc  =  0.15  (8.15) 

The  inflow  and  outflow  boundary  conditions  used  in  this  case  are  described  in 
previous  section.  The  far  field  upper  and  lower  boundaries  are  treated  as  perfect 
non-reflective  outflow  boundary. 

The  flow  field  at  inflow  boundary  is  initially  set  to  be  uniform.  The  direction  of  the 
incoming  flow  is  parallel  to  the  constant  rj  lines.  Obviously,  at  the  inflow  boundary, 
the  transverse  flux  and  viscous  terms  are  very  small.  The  error  caused  by  LODI 
relations  is  very  small  and  negligible. 

To  compare  the  present  NRBC  with  the  commonly-used  outflow  boundary  condi¬ 
tions  (FPBC)  in  which  the  fixed  downstream  pressure  is  imposed  for  subsonic  flow 
simulations,  the  first  computation  is  carried  out  with  FPBC,  which  includes:  at  in¬ 
flow  boundary,  ua,  vrj)  p0  are  given  such  that  the  streamwise  velocity  component  is 
uniformly  distributed  and  the  transverse  velocity  component  is  equal  to  zero,  the 
pressure  is  extrapolated  from  the  interior  domain,  pQ  =  p*.  then  the  total  energy  peQ 
can  be  computed  from  the  equation  of  state;  at  outflow  boundary,  all  the  primitive 
variables  are  extrapolated  from  the  inner  domain  except  the  pressure  is  set  to  be 
constant.  The  CFL  number  of  the  pseudo  time  step  is  500.  The  physical  time  step 
CFL  is  0.74,  which  is  determined  by  the  time  accuracy  of  the  physical  problem 

The  same  case  with  the  same  initial  flow  condition  and  same  CFL  number  is  then 
calculated  using  the  NRBC  developed  in  present  study.  The  value  of  a  in  Equation 
(5.27)  is  set  to  be  0.25.  Figures  12.58  and  12.59  show  the  computed  density  contours  at 
four  instants  by  FPBC  and  NRBC  respectively.  It  can  be  seen  from  figure  12.58  that 
the  flow  field  is  seriously  distorted  by  the  reflective  waves  when  the  vortex  propagates 
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through  the  exit  boundary.  But  there  is  no  noticeable  distortion  in  the  solution 
calculated  using  NRBC  as  shown  in  figure  12.59.  The  vortex  passes  through  the 
outflow  NRBC  very  smoothly.  Figures  12.60  and  12.61  show  the  relative  streamwise 
velocity  component,  (ur  —  uT00)/uT00  contours  at  the  same  four  instants  by  FPBC 
and  NRBC  respectively.  The  same  phenomenon  is  observed.  Figure  12.62  shows  the 
time  histories  of  | pmax  —pmin |  for  the  full  course  of  the  computation.  As  can  be  seen, 
the  level  of  spurious  pressure  reflection  caused  by  FPBC  is  much  higher  than  the  one 
by  NRBC. 

8.5.2  Inlet-Diffuser  Flow 

To  test  the  non-reflective  boundary  conditions  for  realistic  engineering  problems.  A 
transonic  inlet-diffuser  with  shock  wave  boundary  layer  interaction  [62]  is  computed 
to  demonstrate  the  advantage  of  the  NRBC. 


8.5.3  Steady  State  Solutions 

The  steady  state  solution  of  the  inlet  diffuser  is  calculated  first  to  verify  that  the 
NRBC  is  consistent  with  the  steady  state  flow.  The  Reynolds  Number  is  3.45 xlO5 
and  the  inlet  Mach  number  is  0.46. 

The  baseline  geometry  of  the  inlet  diffuser  has  a  height  of  H  =  4.4cm  at  the  throat 
and  a  total  length  of  12 ,6H.  This  case  is  run  using  a  H-type  grid  with  the  dimensions 
of  110x56.  The  turbulence  shear  stress  and  heat  flux  are  calculated  by  the  Baldwin- 
Lomax  model[48].  The  experimental  data  provided  by  Bogar  et  al. [62]  are  available 
for  validation. 

As  discussed  before,  for  3-D  case,  at  £  =  1  inlet  boundary,  waves  L\  -  £4  enter  the 
boundary  and  £5  leaves.  Hence  four  physical  boundary  conditions  are  required  at  this 
boundary.  The  amplitude  of  the  outgoing  characteristic  wave  £5  can  be  estimated 
from  the  interior  points. 

According  to  [33],  the  inlet  and  wall  NRBC  are  not  as  critical  as  the  exit  NRBC.  For 
this  transonic  inlet-diffuser  case,  at  the  upstream  of  the  shock,  the  flow  is  supersonic. 
Hence  the  perturbation  will  not  propagate  upstream.  The  NRBC  at  inlet  therefore 
may  not  be  necessary.  The  conventional  boundary  condition  at  inlet  is  expected  to 
work  well.  However,  at  the  downstream  of  the  shock,  the  flow  is  subsonic.  The 
oscillation  of  the  shock  will  generate  strong  reflecting  waves  at  the  exit  boundary. 
Therefore  exit  NRBC  is  essential  for  this  case.  For  this  reason,  the  inflow  NRBC 
is  not  used  in  this  study.  Instead,  the  inlet  BC  with  given  total  pressure  Pt,  total 
temperature  Tt,  and  flow  angle  is  used.  The  NRBC  outflow  and  wall  conditions  used 
are  those  described  in  previous  section. 

The  Mach  number  contours  are  shown  in  Figure  12.63.  Corresponding  to  different 
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back  pressure  in  the  experiment,  there  are  two  cases  of  the  flow,  one  has  a  weak 
shock  {poUtiet/Pt— 0.82)  and  the  other  has  a  strong  shock  {jpoutiet/Pt= 0.72).  Figures 
12.64  and  12.65  show  the  computed  static  pressure  distributions  compared  with  the 
experimental  data  along  the  top  and  bottom  wall  for  the  weak  shock  case.  Good 
agreement  is  obtained  between  the  computation  and  experiment. 

Figures  12.66  and  12.67  show  the  static  pressure  distribution  compared  with  the 
experimental  data  along  the  top  and  bottom  walls  for  the  strong  shock  case.  Due 
to  the  strong  shock  interacting  with  the  turbulent  boundary  layer,  there  is  a  flow 
separation  downstream  of  the  shock,  which  is  not  well  predicted.  There  may  be  two 
reasons  for  the  problem:  1)  the  flow  is  unsteady  due  to  the  separation  and  hence 
the  steady  state  solution  can  not  capture  the  separation  bubble  length  correctly;  2) 
the  Baldwin-Lomax  turbulence  model  is  inadequate  to  handle  the  non-equilibrium 
separated  flow. 

Different  a  values  from  0.1  to  0.35  are  tested  and  the  results  show  that  the  steady 
state  results  are  insensitive  to  the  a  value.  The  fixed  pressure  boundary  conditions  are 
also  applied  to  the  same  case,  and  achieve  almost  the  same  results  as  those  computed 
by  the  NRBC.  This  is  because  that,  for  the  steady  state  solutions,  the  reflective  waves 
are  eventually  diffused  when  the  steady  state  solution  is  converged. 


8.5.4  Unsteady  Solutions 

The  steady  state  calculation  indicates  that  the  NRBC  is  not  essential  since  the  ar¬ 
tificial  reflective  waves  are  diffused  when  the  solution  is  converged.  However,  it  is 
very  different  when  the  unsteady  flow  is  calculated.  For  the  inlet-diffuser  case  with  a 
strong  shock  wave,  the  fixed  pressure  outflow  boundary  conditions  (FPBC)  generates 
strong  reflective  waves  and  makes  the  shock  wave  severely  oscillating  inside  the  duct. 
The  amplitude  of  the  oscillation  is  far  greater  than  the  experimental  results.  When 
the  NRBC  is  applied,  the  shock  oscillation  is  dramatically  reduced. 

Figures  12.68  and  12.69  show  the  time  averaged  pressure  distributions  compared 
with  the  experimental  data.  The  CFL  number  of  the  pseudo  time  step  is  equal  to  5 
for  all  the  cases.  The  CFL  number  of  the  physical  time  step  is  about  3000,  which  gives 
132  steps  in  a  shock  oscillation  cycle  of  the  computed  dominant  frequency  260(Hz). 
The  value  of  a  in  Equation  (5.27)  is  0.25.  Due  to  the  large  shock  oscillation,  the 
shock  location  is  smeared  out  for  the  fixed  pressure  outflow  boundary  conditions 
(FPBC).  Hence  the  shock  location,  strength  and  the  pressure  downstream  of  the 
shock  are  poorly  predicted.  When  the  NRBC  is  applied,  the  reduced  shock  oscillation 
yields  sharp  shock  profile  and  much  better  agreement  of  the  shock  location  with  the 
experiment. 

To  match  the  experimental  geometry  for  the  measured  shock  oscillation  frequency, 
the  computational  domain  is  then  extend  to  a  total  length  of  21.3 H.  The  simulation 
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is  carried  out  with  FPBC  and  NRBC  respectively  using  the  same  flow  conditions. 
The  computed  and  experimental  power  spectra  for  the  static  pressure  at  the  exit 
location  ( x/H  =  14.218)  are  shown  in  Figure  12.70.  The  experimental  spectrum 
(the  bottom  one)  of  the  shock  oscillation  measured  for  this  case  has  one  dominant 
frequency  around  200  (Hz)[ 62],  The  power  spectrum  computed  by  NRBC  in  the 
middle  plot  shows  that  there  is  only  one  significantly  dominant  frequency  at  about 
260  (Hz).  The  power  spectrum  computed  by  FPBC  in  the  top  plot  has  the  dominant 
frequency  at  380  (Hz),  which  is  obviously  very  different  from  the  dominant  frequency 
of  the  experimental  value  of  200  (Hz).  Obviously,  NRBC  gives  better  results  than 
the  FPBS  does,  which  is  also  better  than  the  computed  value  of  317  (Hz)  predicted 
by  Hsieh  et  al.  [77].  It  is  evident  that  the  NRBC  improves  the  numerical  accuracy 
by  reducing  the  false  reflections,  and  the  noise  level  created  by  NRBC  is  also  much 
lower  than  the  one  created  by  FPBC.  The  results  show  that  the  NRBC  is  essential 
to  accurately  predict  unsteady  aerodynamic  forcing. 


8.6  Separated  Flow  of  NASA  3D  Flutter  Cascade 

NASA  GRC  conducted  compressor  cascade  wind  tunnel  tests  to  mimic  the  flow  field 
of  transonic  compressor  blade  flutter [78]  [79].  The  cascade  is  shown  in  fig.  12.71.  As 
the  first  validation  step,  the  steady  state  solution  of  the  3D  cascade  is  calculated.  The 
3D  effect  of  the  flow  field  is  from  two  aspects:  1)  the  cascade  geometry  is  different 
near  the  end  walls  so  that  the  blades  can  be  mounted;  2)  there  is  no  endwall  boundary 
suction.  Hence  the  endwall  boundary  layer  will  create  the  3D  effect  in  the  spanwise 
direction.  The  parameters  used  in  the  calculation  is  listed  in  Table  8.7. 


Chord 

8.89  cm 

Stagger  angle 

60° 

Solidity 

1.52 

Reynolds  number 

1.522  x  106 

Inlet  Mach  number 

0.5 

Mesh  size 

150  x  60  x  40 

CFL 

5.0 

Table  8.7:  NASA  3D  Flutter  Cascade 

The  3D  mesh  is  shown  in  Figure  12.72.  The  mesh  is  highly  clustered  near  the 
cascade  surface  and  the  end  walls  to  make  sure  y+  is  within  the  range  of  3.  The 
mesh  is  also  clustered  around  leading  edge  and  trailing  edge.  In  Figure  12.73,  the  2D 
mesh  of  the  top,  midspan  and  bottom  planes  are  plotted,  which  indicate  the  geometry 
variation. 

The  boundary  conditions  are  set  as  the  following:  at  subsonic  inlet,  total  pressure 
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and  total  temperature  are  specified;  at  subsonic  outlet,  the  static  pressure  is  specified; 
no  slip  adiabatic  wall  boundary  condition  are  applied  to  bottom,  top  end  walls  and 
the  cascade  surfaces;  periodic  boundary  conditions  are  used  in  the  pitch  direction 
upstream  and  downstream  of  the  cascade  wall  surfaces. 


8.6.1  Steady  state  results 

Though  the  flow  is  separated  and  unsteady  when  the  incidence  is  high,  the  time 
averaged  flow  field  is  calculated  with  the  local  time  stepping  enabled  and  the  dual 
time  stepping  disabled  in  the  code.  The  solution  is  obtained  when  the  calculated 
flowfield  is  unchanged.  The  steady  state  solution  is  also  used  as  the  initial  solution 
for  the  corresponding  unsteady  calculation. 

At  incidence  angle  of  0°,  the  computed  midspan  surface  static  pressure  distribution 
is  in  a  very  good  agreement  with  the  experimental  measurement  [79]  as  shown  in  Figure 
12.74. 

An  incidence  of  10°  is  chosen  for  the  following  numerical  study.  The  van  Leer 
scheme  is  used  to  evaluate  the  inviscid  flux.  Though  the  van  Leer  scheme  is  more 
diffusive  than  the  Roe  schemee,  when  working  with  the  Baldwin-Lomax  turbulence 
model,  it  gives  better  agreement  with  the  experiment  in  the  current  code.  The  Roe 
scheme  predicts  the  separation  larger  than  the  experiment.  The  inlet  Mach  number 
is  obtained  by  adjusting  the  back  pressure  level. 

The  result  of  the  case  with  inlet  Ma= 0.5  is  described  in  this  section  to  demon¬ 
strated  the  characteristics  of  the  3D  separation  flow  field.  More  results  for  Ma— 0.8 
and  Ma— 1.18  can  be  found  in  [80].  In  the  case  of  Ma=0.5,  the  corresponding 
Reynolds  number  is  9.6699xl05.  The  CFL  in  the  Gauss-Seidel  iteration  is  5.0.  The 
calculation  starts  from  a  flowfield  at  rest. 

The  flow  stream  lines  on  the  mid-span  plane  is  shown  in  Fig.  12.75.  The  flow 
exhibits  a  large  separated  region  on  the  suction  surface  that  starts  immediately  at 
the  leading  edge  and  extends  down  to  45%  of  the  blade  chord.  The  flow  pattern  on 
the  suction  surface  is  plotted  on  the  left  in  Fig.  12.76.  The  separation  region  has  a 
parabola  shape,  which  is  approximately  symmetric  about  along  the  blade  mid-span 
line  and  extends  to  the  blade  upstream  corners.  Two  counter  rotating  vortexes  are 
formed  downstream  of  the  blade  leading  edge  corners  on  the  suction  surface  at  its 
two  ends.  The  experiment  visualization  with  dye  oil  technique  is  shown  in  Fig.  12.76 
on  the  right.  The  computation  results  agree  with  the  experiment  fairly  well,  except 
that  the  numerical  results  shows  a  fuller  separation  region  in  span-wise  direction. 

The  mid-span  static  pressure  distribution  is  plotted  and  compared  with  the  exper¬ 
iment  measurement  in  Fig.  12.77.  A  reasonable  agreement  is  achieved.  The  pressure 
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is  expressed  as  the  pressure  coefficient, 


2  PinUin 

where  p  is  the  local  static  pressure.  pin ,  and  Uin  are  the  averaged  static  pressure, 
density  and  velocity  at  the  inlet. 

On  the  leading  edge  of  the  suction  side,  the  numerical  pressure  results  varies  more 
steeply  than  that  given  by  the  measurement.  The  separation  region  denoted  by  the 
cross  of  the  pressure  distributions  on  suction  and  pressure  surfaces  agrees  very  well 
with  the  experiment. 

The  separation  region  is  enlarged  when  Roe  scheme  is  used  to  calculate  the  inviscid 
flux.  A  possible  reason  for  the  difference  is  the  application  of  turbulence  model  on 
the  H-type  mesh  in  the  current  code.  It  is  shown  that  the  implementation  details  the 
Baldwin-Lomax  model  is  vital  to  the  resulted  turbulent  viscosity  accuracy.  Better 
agreement  is  obtained  in  [81]  for  the  same  case,  where  the  Roe  scheme  is  applied 
on  an  O-type  mesh.  The  H-type  mesh  in  the  current  code  is  generated  by  an  ellip¬ 
tic  method  as  a  whole,  which  meets  difficulty  in  the  mesh  orthogonality  in  the  wall 
boundary  region,  which  affects  the  accuracy  in  calculating  the  outer  layer  eddy  vis¬ 
cosity  coefficient.  A  two-layer  H-type  mesh  is  used  in  [82],  where  an  inner  algebraic 
mesh  is  surrounded  by  an  outer  elliptic  method  generated  mesh.  The  the  inner  mesh 
is  designed  to  achieve  better  orthogonality.  These  grid  generation  techniques  will  be 
implemented  in  the  code  in  the  next  step  research  work. 

8.6.2  Unsteady  separated  flow  simulation 

The  separation  is  believed  to  bring  high  unsteadiness  to  the  cascade  flow  pattern.  The 
inlet  Mach  number  is  an  important  factor  which  affects  the  separation  characteristics  [80] 
To  study  the  influence  of  the  inlet  Mach  number  on  the  unsteady  characteristics  of 
the  separated  flow,  numerical  simulation  is  carried  out  for  high  incidence  angle  cases 
with  Mach  number  0.5,  0.8  and  1.18.  Each  unsteady  calculation  is  carried  out  based 
on  its  corresponding  steady  state  result. 

Due  to  the  limitation  of  the  computation  capability,  the  physical  time  interval  is 
chosen  as  large  as  10%  of  the  characteristic  time  of  the  cascade,  tc  =  c/Uin.  f7j„is  the 
inlet  velocity.  This  time  interval  varies  with  the  inlet  Mach  number.  The  CFL  number 
used  in  the  pseudo  time  Gauss-Seidel  implicit  iteration  is  20.0.  Twenty  pseudo  time 
steps  are  used  for  each  physical  time  step.  Two  parameters  are  recorded  to  analyze  the 
unsteady  characteristics  of  the  separation  flow.  The  first  is  the  mid-span  separation 
bubble  length(x),  which  is  marked  by  the  streamwise  zero  velocity  point  at  the  first 
inner  mesh  point  on  the  suction  surface.  The  second  parameter  is  the  unsteady  static 
pressure  (p)  measured  at  the  location  of  13%  downstream  the  leading  edge  on  the 
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suction  surface,  which  is  the  same  as  the  experiment  measurement  location.  The 
unsteady  pressure  is  refereed  as  “check  point  pressure”  in  the  following. 


Ma  =  0.5 

In  the  case  of  Ma  =  0.5,  the  physical  time  interval  is  set  as  0.052744  ms. 

Fig.  12.78  shows  the  time  history  of  the  separation  length  oscillation  in  a  time 
segment  of  16  ms  (30fc).  The  separation  length  increases  rapidly  from  0.45c  to  0.66c 
in  the  first  1.63  ms  (3fc)  and  then  decreases  to  0.63c  at  t  =2  ms  (3.8tc).  The  separation 
region  then  grows  up  again  toward  downstream  to  0.73c  at  t  =3  ms  (5.7 tc).  With 
the  time  progressing,  the  separation  region  boundary  oscillates  back  and  forth  on  the 
suction  surface.  The  average  length  tends  to  increase  gradually  until  a  periodic  state 
is  reached  after  4.4  ms  (8.3fc).  The  oscillation  of  the  separation  length  is  between 
0.73c  and  0.76c  with  a  fixed  cycle.  The  periodicity  information  is  clearly  extracted 
using  the  FFT  technique.  The  separation  region  oscillation  spectra  is  calculated  from 
the  unsteady  data  after  4.4  ms  (8.3fc).  The  frequency  spectrum  is  shown  in  Fig.  12.79, 
which  clearly  shows  a  peak  at  770  Hz.  This  indicates  the  separation  length  oscillates 
with  a  period  of  1.25  ms  (2.37fc). 

Compared  with  the  steady  state  solution,  the  unsteady  separation  calculation  re¬ 
sults  in  larger  separation  size.  The  reason  is  not  clear. 

The  unsteady  check  point  pressure  data  shows  similar  characteristics  of  the  un¬ 
steady  separation  flow.  Fig.  12.80  shows  a  segment  of  16  ms  (30tc)  pressure  oscillation 
data.  The  pressure  start  at  p  =2.659  from  the  steady  state  results.  The  oscillation 
is  between  p  =2.68  and  p  =2.72  after  t  =4.4  ms  (8.3tc).  The  oscillation  amplitude 
is  about  1.5%  of  the  averaged  pressure  level.  The  frequency  spectrum  is  shown  in 
Fig.  12.81  with  a  peak  at  770  Hz. 

The  mechanism  behind  the  unsteady  characteristics  of  the  separation  is  illustrated 
in  Fig.  12.83,  where  the  evolution  of  a  separation  oscillation  cycle  is  plotted.  The 
stream  lines  at  8  time  steps  show  the  leading  edge  vortex  shedding  development. 
There  are  4  physical  time  steps  (0.21  ms,  0.4ic)  between  2  sequential  plots.  The 
relationship  between  the  oscillation  of  the  pressure  and  the  separation  length  is  shown 
in  Fig.  12.82. 

At  the  starting  time  level  a  ( t  =4.4305  ms,  8.4fc),  the  separation  region  has  just 
passed  the  maximum  length  location.  There  are  2  vortexes  in  the  separation  bubble. 
They  are  rotating  in  the  same  direction.  The  check  point  pressure  is  going  up.  At 
time  level  b  ( t  =4.6415  ms,  8.8fc),  the  two  vortexes  are  pushed  toward  downstream. 
The  second  vortex  diminishes.  The  first  vortex  grows  quickly  and  becomes  the  only 
vortex  in  the  separation  bubble.  The  separation  region  becomes  thicker  in  the  span- 
wise  direction,  but  shorter  in  the  stream-wise  direction.  The  surface  check  point 
pressure  reaches  its  maximum  level  at  this  time  level.  At  time  level  c  ( t  =4.8524  ms, 
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9.2fc),  the  separation  bubble  approaches  its  shortest  length  in  stream- wise  direction, 
and  maximum  thickness  in  the  span-wise  direction.  The  check  point  pressure  is 
going  down.  At  time  level  d  (t  —  5.0634ms,  9.6tc),  the  separation  bubble  has  passed 
the  minimum  length  location,  and  begins  to  extend  toward  downstream.  The  check 
point  pressure  is  still  going  down.  At  time  level  e  ( t  =5.2744  ms,  10.0 tc),  a  new 
vortex  is  generated  at  leading  edge  and  becomes  the  first  vortex.  The  check  point 
pressure  reaches  its  minimum  value.  The  separation  length  is  still  increasing.  At  time 
level  /  (t=4854  ms,  10.4tc),  the  first  vortex  continues  to  grow.  The  second  vortex 
is  pushed  toward  downstream.  Both  the  check  point  pressure  and  the  separation 
length  are  going  up.  The  latter  is  approaching  its  maximum  location.  At  time  level  g 
(t  =5.6964  ms,  10.8tc),  the  two  vortexes  have  almost  the  same  size,  the  flow  structure 
is  close  to  the  starting  time  level  a.  The  separation  boundary  has  passed  its  maximum 
location  and  begin  to  shrink  toward  upstream.  The  check  point  pressure  is  going  up. 
At  time  level  h  ( t  =5.9073  ms,  11.2tc),  the  second  vortex  diminishes.  A  new  cycle  is 
started  at  this  time  level. 

The  leading  edge  vortex  shedding  exhibits  obvious  periodical  pattern  in  its  evolu¬ 
tion  process.  The  leading  edge  keeps  generating  new  vortexes.  The  new  vortex  pushes 
the  old  vortex  bubble  toward  downstream  and  the  old  vortex  decreases  in  size  at  the 
same  time.  When  the  two  vortexes  become  of  the  same  size,  the  maximum  separa¬ 
tion  length  is  reached,  where  the  separation  bubble  has  the  thinnest  size  in  span-wise 
direction.  As  the  new  vortex  grows  further,  the  old  vortex  will  diminish.  The  sepa¬ 
ration  bubble  will  move  upstream  and  makes  the  bubble  thicker.  The  leading  edge 
surface  check  point  pressure  reaches  its  maximum  level  when  the  separation  bubble 
shrinks  and  reaches  its  minimum  level  when  the  separation  bubble  boundary  extends. 
The  vortex  generation,  pressure  variation  and  separation  length  oscillation  have  the 
same  frequency  characteristics  with  a  phase  difference  as  shown  in  Fig.  12.82.  Such 
oscillation  is  maybe  one  of  the  reasons  to  cause  to  flutter. 


Ma  =  0.8 

The  physical  time  interval  used  in  the  calculation  for  the  case  of  Mach  number  0.8  is 
0.0342  ms.  A  similar  flowfield  unsteady  characteristics  is  exhibited  in  the  computation 
results. 

Fig.  12.84  and  Fig.  12.86  show  the  separation  length  and  the  checkpoint  static 
pressure  oscillation  history  in  a  time  period  of  21  ms.  A  clear  periodicity  is  shown 
in  these  two  figures.  It  is  found  in  the  time  averaged  steady  state  study  in  [80]  that, 
the  increase  of  the  inlet  Mach  number  will  enlarge  the  separation  bubble  in  size. 
Fig.  12.86  indicates  that  the  inlet  Mach  number  increase  also  increases  the  amplitude 
of  the  pressure  oscillation.  The  oscillation  amplitude  is  increased  to  about  5%  of 
the  averaged  pressure  level.  The  increased  kinetic  energy  in  the  inflow  bring  higher 
unsteadiness  intensity  to  the  separated  flow  filed. 
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The  corresponding  FFT  frequency  analysis  is  shown  in  Fig.  12.85  and  Fig.  12.87 
respectively.  The  frequency  analysis  is  based  on  the  oscillation  data  after  f=5  ms. 
The  unsteady  separation  flow  exhibits  higher  oscillation  frequency  because  of  the 
increased  inlet  Mach  number.  A  clear  frequency  spectrum  peak  is  shown  at  1400  Hz 
in  both  figures,  twice  the  frequency  in  the  case  of  Ma  =  0.5. 


Ma  =  1.18 

In  the  steady  state  simulation  of  the  cascade  at  Ma= 1.18  in  [80],  the  separation 
flow  characteristics  are  very  different  from  those  at  subsonic.  The  further  increased 
kinetic  energy  in  the  inflow  makes  the  flow  attached  to  the  blade  surface  in  the 
leading  edge.  A  smaller  sized  separation  region  appears  after  the  shock  wave  because 
of  the  interaction  of  the  shock  wave  and  the  turbulent  boundary  layer.  The  separation 
bubble  shrinks  in  size  and  is  located  only  in  a  small  region  at  the  center  of  the  suction 
surface  region. 

The  physical  time  interval  in  the  calculation  is  set  as  0.02483  ms.  The  pressure 
check  point  is  located  outside  of  the  separation  region  in  the  supersonic  case.  The 
pressure  oscillation  history  is  shown  in  Fig.  12.88.  The  oscillation  amplitude  is  very 
small  compared  to  the  cases  of  Ma  =  0.5  and  Ma  =  0.8.  The  flow  tends  to  steady  at 
the  leading  edge.  The  pressure  oscillation  frequency  spectrum  is  shown  in  Fig.  12.89. 

The  computed  characteristics  of  the  separation  flow  above  is  similar  to  the  experi¬ 
ment  measurement  in  [83].  In  [83],  when  the  blade  is  fixed,  the  blade  surface  pressure 
for  low  subsonic  inlet  flow  at  Ma= 0.5  and  low  supersonic  inlet  flow  at  Ma= 1.1  ex¬ 
hibits  very  low  unsteadiness  and  very  strong  self-induced  oscillations  with  a  frequency 
of  110Hz  is  observed  in  high  subsonic  inlet  flow  at  Ma= 0.8.  However,  the  strong  low 
frequency  oscillation  is  attributed  to  the  tunnel  resonance  characteristics  instead  of 
the  flow  unsteadiness  due  to  the  flow  separation  in  [84].  Even  though,  the  cascade 
flow  separation  is  believed  to  have  a  direct  relation  with  the  wall  surface  pressure 
unsteady  oscillation,  which  is  an  important  factor  to  the  flutter.  Further  detailed 
numerical  research  is  necessary  to  discover  the  mechanism. 


8.7  Forced  Vibration  of  the  NASA  Flutter  Cas¬ 
cade 


8.7.1  Parameters  used  in  flow  analysis 


The  steady  state  static  pressure  coefficient  is  defined  as  the  following, 


Cp  (x/c) 


p{x/c)  —  Pop 
\PooUl 


(8.16) 
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where  p^,  p^  and  Uoo  are  the  averaged  pressure,  density  and  velocity  at  inlet. 

The  blades  vibrate  harmonically  with  a  constant  IBPA.  The  motion  of  the  nth 
blade  is  defined  by  the  blade  deflection  angle  [?], 

an  ( t )  =  c*o  +  &Re  [exp  ( i  (ut  +  n{3 ))]  (8-17) 


where  n  is  the  blade  index,  t  is  the  time,  ao  is  the  deflection  angle  at  the  mean  blade 
position,  a  is  the  amplitude  of  blade  deflection,  Re  denotes  the  real  part  of  a  complex 
value,  u)  is  the  angular  frequency,  (3  is  the  inter  blade  phase  angle. 

The  reduced  frequency  kc  is  defined  based  on  the  chord  length  C  as  the  following, 


kc  — 


LoC 

un a 


(8.18) 


The  first  harmonic  unsteady  pressure  coefficient  is  defined  as, 


Cv  (x)  = 


Pi  (s) 

PooU^a 


(8.19) 


where  px  (x)  is  first  harmonic  pressure  along  the  blade  surface.  It  is  a  complex  value 
obtained  from  the  unsteady  pressure  signals  using  the  Fourier  transformation,  pi  ( x ) 
has  a  phase  angle  relative  to  the  blade  motion  a. 

The  time  dependent  aerodynamic  moment  coefficient  is  defined  as, 


Cm(t)  = 


(8.20) 


where  p  ( x )  is  the  unsteady  pressure  along  the  blade  surface,  s  is  the  surface  area  vector 
pointing  outward  from  the  blade,  r  is  the  vector  pointing  from  the  pivot  location  to 
an  arbitrary  point  x  on  the  surface. 

The  imaginary  part,  or  the  out  of  phase  part  of  the  unsteady  pressure  determines 
the  damping  or  excitation  of  the  blade  motion.  The  aerodynamic  damping  coefficient 
is  defined  as, 

E  =  -Im(Cm)  (8.21) 

where  Im  denotes  the  imaginary  part  of  a  complex  value.  A  positive  E  corresponds 
to  a  damped  oscillation. 


8.7.2  The  Cascade 

The  NASA  Lewis  Oscillating  Cascade  test  section  consists  of  9  identical  airfoils  with 
the  cross  section  similar  to  the  tip  airfoil  of  modern  low- aspect  ratio  fan  blades  [?]. 
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The  airfoil  has  a  chord  of  8.89  cm  and  is  installed  with  a  stagger  angle  of  60°.  The 
solidity  is  1.52.  In  the  experiment,  the  inlet  Mach  number  is  0.5.  All  blades  vibrate 
simultaneously  along  a  pitching  axis  at  0.5  chord  with  a  constant  IBPA  of  180°.  The 
oscillating  amplitude  is  1.2°  and  the  reduced  frequency  based  on  chord  varies  at  0.4, 
0.8  and  1.2.  The  blade  motion  is  identical  on  every  other  blade.  The  two  neighboring 
blades  always  vibrate  in  opposite  direction. 

The  inlet  Mach  number  Ma  =  0.5  is  achieved  by  adjusting  the  outlet  static  pressure 
level.  The  Reynolds  number  based  on  the  chord  length  is  9x  105.  The  flow  incidence  is 
0°.  For  the  unsteady  dual- time  stepping,  one  physical  blade  oscillation  cycle  is  divided 
into  100  time  intervals  and  100  pseudo  time  Gauss-Seidel  iterations  are  carried  out  for 
each  physical  time  step.  The  100  pseudo  time  iterations  are  proved  to  be  sufficient  to 
obtain  a  converged  solution  within  a  physical  time  step  with  the  residual  reduced  by  3 
orders  in  magnitude.  Before  the  unsteady  simulation,  the  corresponding  steady  state 
calculation  is  carried  out  to  obtain  the  initial  flow  field  for  the  unsteady  computation. 
The  inlet  flow  angle  is  adjusted  to  61°  to  obtain  good  agreement  with  the  steady  state 
experimental  surface  pressure  distribution. 

Before  the  full  scale  computation,  a  simplified  2-passage  cascade  is  computed  with 
a  reduced  frequency  of  kc  =  0.8  by  applying  the  periodic  boundary  condition.  The 
information  of  the  moving  mesh  and  the  conservative  variables  are  exchanged  across 
the  periodic  boundary.  The  multi-passage  full  scale  simulation  is  then  conducted 
for  the  9-blade  cascade  with  the  wind  tunnel  side  walls  for  more  realistic  results. 
The  periodic  boundary  condition  is  not  needed.  The  vibration  frequencies  of  fcc=0.4, 
0.8  and  1.2  are  calculated  and  compared  with  the  experiment  measurement.  The 
influence  of  the  end  walls  are  studied  by  comparing  the  results  of  A;c=0.8  from  both 
simulations. 


8.7.3  Computation  domain  decomposition  and  mesh  gener¬ 
ation 

As  shown  in  Fig.  12.107,  the  computation  domain  which  is  consistent  with  the  exper¬ 
iment  configuration  is  split  into  10  subdomains  based  on  flow  passages  PI  through 
P10.  The  10  subdomains  are  computed  by  10  CPUs  running  in  parallel.  The  subdo¬ 
mains  are  separated  from  their  neighbors  by  the  blade  surfaces  B1  through  B9  and  the 
MPI  interface  boundaries.  The  MPI  interface  boundaries  are  straight  lines  passing 
through  the  leading  edge  (LE)  and  trailing  edge  (TE)  of  the  blades  with  appropriate 
angles  in  accord  with  the  local  flow  direction.  The  pitchwise  distance  between  the 
end  wall  and  the  blade  is  half  of  the  inner  pitch  distance.  The  US  (upper  surface)  is 
the  suction  surface  and  LS  (lower  surface)  is  the  pressure  surface. 

The  inlet  and  outlet  boundaries  are  set  as  1.5  and  3  times  chord  length  away  from 
the  airfoil  LE  and  TE  in  axial  direction.  A  part  of  the  mesh  is  shown  in  Fig.  12.108 
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with  the  regions  of  LE  and  TE  zoomed  in  for  more  details.  The  H-type  mesh  is 
generated  for  each  subdomain  respectively.  Each  subdomain  shares  the  grid  point 
distribution  on  the  common  MPI  interface  boundaries  with  its  neighbors.  To  achieve 
good  orthogonality  on  the  blade  surface,  an  additional  algebraic  boundary  layer  mesh 
is  generated  in  the  wall  surface  region.  As  shown  in  Fig.  12.108,  the  grid  lines  are 
orthogonal  on  all  blade  surfaces  except  the  small  regions  at  LE  and  TE.  The  mesh 
size  is  195(£)  x  180(7?)  for  all  subdomains.  For  clarity,  the  mesh  is  plotted  every  4  lines 
in  the  un-zoomed  plots.  The  blade  surface  has  100  points  in  streamwise  direction. 
The  boundary  layer  has  40  points  in  pitchwise  direction.  Because  of  the  high  gradient 
of  the  flow  variables  in  near  wall  region,  the  mesh  is  clustered  near  the  wall  surfaces. 
On  the  blade  surfaces,  the  grid  points  are  also  clustered  toward  the  LE  and  the  TE 
in  streamwise  direction. 

8.7.4  Simulation  in  two  passage  cascade 

The  two  passage  cascade  simulation  uses  the  meshes  of  two  inner  neighboring  passages 
(P2  and  P3)  in  the  compressor  cascade.  As  shown  in  Fig.  12.109,  The  blade  between 
the  two  passages  is  called  BC  (blade  at  center)  and  the  two  blade  surfaces  on  the  two 
outside  periodic  boundaries  are  treated  as  blade  BP  (blade  at  periodic  boundary)  . 

The  steady  state  pressure  coefficient  distributions  along  the  blade  surfaces  are 
plotted  in  Fig.  12.110.  The  pressure  coefficient  predicted  agrees  well  with  the  mea¬ 
surement.  The  result  on  BC  agrees  very  well  with  that  of  BP,  which  shows  good 
periodicity  is  achieved  on  the  pitchwise  direction. 

Fig.  12.111  shows  the  pressure  variation  history  on  two  points  on  the  suction 
surface  (US)  and  the  pressure  surface  (LS)  respectively.  The  pressure  on  the  suction 
surface  is  located  at  x/C  =  0.15  and  the  pressure  on  the  pressure  surface  is  located 
at  x/C  =  0.1.  The  temporal  periodicity  is  achieved  very  soon  after  the  start  of  the 
vibration  simulation.  Because  of  the  excellent  temporal  periodicity,  the  unsteady 
data  extracted  from  a  single  blade  motion  cycle  is  enough  for  the  unsteady  Fourier 
analysis.  The  IBPA  of  180°  is  clearly  shown  by  comparing  the  pressure  maximums 
and  minimums  on  BC  and  BP. 

The  unsteady  pressure  coefficients  are  plotted  in  Fig.  12.112  and  Fig.  12.113  for 
suction  surface  (US)  and  the  pressure  surface  (LS)  respectively.  The  unsteady  pres¬ 
sure  coefficient  Cp  is  expressed  in  terms  of  the  real  part  or  in  phase  part  and  the 
imaginary  part  or  out  of  phase  part.  On  the  suction  surface,  as  shown  in  Fig.  12.112, 
the  CFD  results  compare  fairly  well  with  experiment  data  after  30%  chord.  The  real 
part  of  the  coefficient  is  predicted  lower  than  experiment  data  on  leading  edge,  but 
the  trend  agrees  very  well  with  the  experiment.  The  imaginary  part  is  over  predicted 
in  the  leading  edge  region.  On  the  pressure  surface,  as  shown  in  Fig.  12.113,  the 
real  part  of  the  unsteady  pressure  coefficient  agrees  well  with  the  experiment  data. 
The  imaginary  part  is  under  predicted  compared  with  the  measurement  on  the  front 
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part  of  the  blade.  This  means  that  the  CFD  does  not  accurately  capture  the  phase 
angle  difference  between  the  pressure  response  and  the  blade  motion.  A  local  flutter 
stability  analysis  based  on  the  aerodynamic  work  per  cycle  suggested  by  Buffum[?] 
is  presented  in  Fig.  12.114  by  plotting  (0.5  -  Im(CPtU pper  —  CPt\ ower)lst-  The  cur¬ 
rent  numerical  simulation  predicts  a  larger  local  stable  region  on  the  front  part  of 
the  blade.  On  the  aft  part,  the  experiment  data  indicates  a  shallow  stable  region. 
The  CFD  predicts  the  trend  very  well.  The  aft  part  stable  region  is  predicted  more 
shallow  than  the  experiment. 


8.7.5  Simulation  in  full  scale  cascade 

The  full  scale  steady  and  unsteady  simulation  use  the  same  inlet  flow  angle  as  that 
of  the  2-passage  case.  In  the  steady  state  calculation,  all  blade  are  parallel  to  each 
other  at  their  mean  positions.  As  shown  in  Fig.  12.107,  the  end  wall  is  made  up  of  3 
sections,  which  have  different  angles  relative  to  the  x  axis:  £*i  =  61°,  «2  =  60°  and 
c*3  =  64°.  The  middle  section  is  parallel  to  the  blade  at  its  mean  position.  The  front 
and  aft  sections  follow  the  inlet  and  outlet  averaged  flow  directions  obtained  in  the 
2-passage  steady  state  computation. 

The  steady  state  Mach  number  contours  for  the  full  scale  cascade  is  shown  in 
Fig.  12.115.  The  flow  pattern  is  highly  influenced  by  the  end  wall  especially  for 
the  near  end  wall  passages,  PI  and  P10.  The  influence  is  reduced  rapidly  from 
the  boundary  passages  to  the  inner  passages.  Good  periodicity  in  flow  pattern  is 
achieved  among  the  inner  passages  (P3  through  P8).  Three  center  blades,  B4,  B5 
and  B6  are  chosen  to  study  the  steady  and  unsteady  periodicity  in  the  rest  of  the 
paper.  Even  though  the  periodicity  looks  good  in  the  Mach  contour  plot,  the  static 
pressure  distribution  still  shows  the  influence  of  the  end  walls  on  different  blades. 

Fig.  12.116  shows  the  pressure  coefficient  chordwise  distribution  on  the  3  center 
blades.  The  experiment  measurement  and  the  2-passage  calculation  results  are  also 
plotted  for  comparison.  The  surface  pressure  increases  gradually  from  B4  to  B6  on 
both  the  pressure  surface  and  the  suction  surface.  The  pressure  distribution  on  blade 
B6  is  closest  to  the  2-passage  periodic  results  on  most  part  of  the  surfaces.  The 
experiment  measurement  also  shows  the  pressure  variation  on  different  bladesf?],  but 
its  variation  trend  is  opposite  to  the  current  numerical  results.  A  possible  reason  for 
this  difference  is  that  the  inlet  and  outlet  end  wall  angles  used  in  the  experiment  may 
differ  from  the  values  used  in  the  current  simulation.  The  experimental  angles  are  not 
available.  Such  a  pitchwise  flow  pattern  difference  is  also  expected  in  the  following 
unsteady  calculations. 

The  full  scale  unsteady  simulation  is  first  carried  out  for  a  reduced  frequency  kc 
=  0.8  to  study  the  end  wall  influence  on  the  periodicity  of  the  blade  unsteady  char¬ 
acteristics.  Figs.  12.117  and  12.118  are  the  unsteady  pressure  coefficient  chordwise 
distribution  on  the  3  center  blades  compared  with  the  2-passage  cascade  and  the  ex- 
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periment  results.  On  the  upper  surface,  as  shown  in  Fig.  12.117,  the  3  blades  have 
very  similar  unsteady  coefficients  on  most  of  the  chordwise  distance.  The  results  of 
blade  B6  are  closest  to  those  of  the  2-passage  calculation.  The  difference  between  the 
full  scale  results  and  the  2-passage  results  mainly  locate  at  the  front  and  center  part 
of  the  blade.  The  full  scale  results  agree  with  the  experiment  better  in  the  center 
part.  On  the  lower  surface,  as  shown  in  Fig.  12.118,  the  full  scale  results  of  the  3 
blades  are  similar.  The  2-passage  results  are  closer  to  the  experiment  data  in  the  real 
part. 

As  shown  in  Fig.  12.119,  the  full  scale  calculations  predict  higher  stability  on  the 
front  part  of  the  blade  compared  with  the  2-passage  results.  The  full  scale  results  are 
closer  to  the  experiment  measurement  on  the  aft  part  of  the  blade.  In  the  chordwise 
region  of  x/C  =  0.5  to  x/C= 0.7,  the  measured  stability  is  better  predicted  in  the  full 
scale  results.  The  2-passage  results  shows  instability  in  the  same  region.  The  end  wall 
influence  on  the  flow  pattern  periodicity  is  clearly  shown  in  the  unsteady  aerodynamic 
moment  oscillation  plots  within  a  whole  blade  motion  cycle  in  Fig.  12.120.  The 
moment  is  plotted  versus  the  normalized  deflection  angle,  a'  =  [a  —  oto )  /ot.  Because 
of  the  end  wall  influence,  the  moment  oscillations  on  the  3  center  blades  are  different. 
They  are  also  different  from  the  2-passage  calculation  results.  The  anti-clockwise 
direction  of  all  the  unsteady  moment  curves  indicates  negative  work  acted  by  the  fluid 
on  the  blade.  The  blade  motion  is  therefore  damped  down  by  the  fluid  flow.  The 
blade  motion  is  stable,  which  corresponds  to  a  positive  damping  coefficient  E.  The 
area  enclosed  by  the  moment  curve  indicates  the  magnitude  of  the  work  exchanged 
between  the  fluid  and  the  blade,  which  is  also  proportional  to  the  magnitude  of  the 
damping  coefficient. 

The  damping  coefficients  on  all  the  9  blades  in  the  full  scale  calculation  versus 
the  blade  number  are  shown  in  Fig.  12.121.  The  damping  coefficient  varies  among 
the  blades.  The  damping  coefficients  for  blade  B4,  B5  and  B6  are  0.67,  0.65  and 
0.68  respectively.  Blade  B1  has  the  lowest  stability  (S  =0.45)  and  blade  B9  has  the 
highest  stability  (E  =  1.4).  The  damping  coefficient  distribution  is  more  uniform  on 
the  center  blades  (B3  through  B7),  even  though  small  variation  exists.  The  blade 
stability  in  the  full  scale  cascade  depends  on  the  location  of  the  blade.  The  damping 
coefficient  obtained  in  the  2-passage  periodic  computation  is  0.55. 

Fig.  12.122  plots  a  series  of  Mach  number  contours  around  blade  B5  and  B6.  A 
separation  bubble  is  generated  and  grows  periodically  on  the  leading  edge  of  the  suc¬ 
tion  surface.  At  t  =0,  the  two  blades  are  initially  located  at  their  mean  positions 
and  are  parallel  to  each  other.  Blade  B5  then  rotates  in  the  counter-clockwise  di¬ 
rection  with  a  negative  deflection  angle  (nose  down).  At  the  same  time,  blade  B6 
is  rotating  in  the  clockwise  direction  with  a  positive  deflection  angle  (nose  up).  At 
f=0.2T,  blade  B5  is  close  to  its  minimum  deflection  position.  The  separation  bubble 
at  its  LE  is  pushed  downstream  and  shrinks  in  size.  At  f=0.4T,  blade  B5  is  rotating 
back  from  its  minimum  deflection  location  toward  its  mean  position,  the  separation 
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bubble  disappears  from  the  suction  surface.  At  t  =  0.8T,  a  new  separation  bubble 
is  generated  when  blade  B5  passes  its  maximum  deflection  position  and  rotates  back 
toward  its  mean  position.  The  bubble  obtains  its  maximum  size  when  the  blade  is 
close  to  its  mean  position.  Similar  phenomenon  is  observed  on  the  neighboring  blade 
B6,  but  with  a  phase  difference  of  180°. 

More  extensive  unsteady  simulations  are  carried  out  for  reduced  frequencies  kc  = 
0.4  and  kc  =  1.2.  Figs.  12.123  and  12.124  show  the  unsteady  pressure  coefficient 
chordwise  distribution  of  kc  —  1.2.  Similar  to  the  results  of  kc  =  0.8,  the  predicted 
unsteady  complex  pressure  coefficients  are  close  to  each  other  on  the  3  center  blades, 
even  though  some  small  difference  exists.  As  shown  in  Fig.  12.123,  on  the  upper 
surface,  the  imaginary  parts  of  CFD  results  agree  very  well  with  the  experiment 
results  except  that  it  is  over-predicted  in  the  region  of  x/C=0.15  to  x/C=0A0.  The 
real  part  is  also  predicted  quite  well  on  the  middle  and  aft  part  of  the  blade.  On 
the  lower  surface,  as  shown  in  Fig.  12.124,  the  predicted  real  parts  compare  very  well 
with  experiment.  The  Imaginary  part  is  under-predicted  from  LE  to  x/C  =  0.7. 

The  local  stability  analysis  for  kc  —  1.2  is  plotted  in  Fig.  12.125.  The  correct 
trend  is  predicted  compared  with  the  experiment  measurement,  even  though  the 
magnitude  does  not  agree  very  well.  The  stability  is  over-predicted  in  LE  region. 
The  unstable  region  predicted  on  the  front  part  of  the  blade  is  smaller  than  the 
experiment  results.  The  stability  is  predicted  on  the  aft  part,  but  the  magnitude  is 
smaller  than  the  experiment  data.  The  damping  coefficients  for  all  the  blades  are 
plotted  in  Fig.  12.126.  The  damping  coefficients  for  blade  B4,  B5  and  B6  are  0.81, 
0.78  and  0.84  respectively.  The  stability  increases  with  the  frequency.  Similar  to  the 
results  of  kc  —  0.8,  the  most  stable  blade  is  blade  B9  (5=1.5)  and  the  least  stable 
blade  is  blade  B1  (E=0.6).  The  damping  coefficient  is  more  uniformly  distributed 
on  the  central  blades.  The  variation  of  the  damping  coefficient  on  the  center  blades 
increases  with  the  increasing  frequency. 

Because  of  the  lack  of  experiment  data,  the  unsteady  pressure  coefficient  of  A;c=0.4 
is  not  presented  for  comparison.  However,  the  local  stability  is  analyzed  and  com¬ 
pared  with  the  experiment  data  in  Fig.  12.127.  Similar  to  the  results  of  kc  =  0.8  and 
kc  =  1.2,  the  trend  is  predicted  well,  but  the  magnitude  differs  from  the  experiment. 
As  expected,  the  damping  coefficient  distribution  is  more  uniform  on  center  blades 
(Fig.  12.128).  The  variation  of  their  magnitudes  decreases  with  the  decreasing  vibra¬ 
tion  frequency  compared  with  the  high  frequency  cases  of  kc  =  0.8  and  kc  =  1.2.  The 
damping  coefficients  on  blade  B4,  B5  and  B6  are  0.448,  0.446  and  0.434  respectively. 
The  most  stable  blade  is  B9  with  5  =  1.02  and  the  least  stable  blade  is  B1  with  5  = 
0.28. 

The  unsteady  aerodynamic  moment  oscillations  on  blade  B5  under  the  3  frequen¬ 
cies  under  investigation  are  plotted  together  and  compared  in  Fig.  12.129.  The  damp¬ 
ing  coefficient  increase  with  the  increasing  frequency  is  indicated  by  the  increased 
area  enclosed  by  the  unsteady  moment  oscillation  curve.  The  local  stability  analysis 
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is  summarized  for  all  the  3  frequencies  in  Fig.  12.130.  The  computation  results  in¬ 
dicate  higher  stability  near  the  leading  edge  for  higher  frequency  vibration,  which  is 
consistent  with  the  experiment  measurement.  Even  though  the  destabilization  region 
on  the  front  part  of  the  blade  and  the  stability  magnitude  on  the  aft  part  of  the  blade 
predicted  by  the  numerical  computation  are  smaller  than  those  in  the  experiment, 
the  trend  is  predicted  well.  Both  the  destabilization  and  stabilization  increase  with 
the  increasing  frequency. 


8.8  3D  AGARD  Wing  Flutter  Prediction 

The  result  of  steady  state  transonic  ONERA  M6  wing  is  calculated  first  in  order 
to  validate  CFD  solver.  Then,  the  flutter  boundary  of  an  AGARD  wing  445.6  is 
calculated. 


8.8.1  Steady  State  Transonic  ONERA  M6  wing 

As  a  validation  of  the  three  dimensional  solver  for  a  transonic  wing,  the  steady  state 
solution  of  the  transonic  ONERA  M6  wing  is  calculated.  The  freestream  conditions 
for  this  study  are  listed  in  Table  8.8  below. 

Table  8.8:  Free-stream  condition  for  ONERA  M6  wing 


Mach  number 

0.8395 

Static  Pressure  (psia) 

12.2913 

Temperature  (R) 

447.0 

Angle-of-Attack  (deg) 

0.0 

Reynolds  Number 

19.7xl06 

This  case  is  calculated  using  an  O-type  grid  with  the  dimension  of  144  (around 
wing)  x60  (normal  to  the  wing)  x40  (spanwise).  The  far  field  boundary  is  located 
15  chords  away  from  the  chord  center  of  the  wing.  The  surface  mesh  of  the  wing  is 
depicted  in  Figure  12.133. 

The  Zha  E-CUSP2  scheme  is  used  to  evaluate  the  inviscid  fluxes  with  the  3rd 
order  MUSCL  type  differencing  [85].  The  turbulent  Reynolds  stress  and  heat  flux  are 
calculated  by  the  Baldwin-Lomax  algebraic  model  [48]. 

The  computed  surface  pressure  distributions  at  various  cross  sections  are  shown  in 
Figure  12.134,  together  with  the  experimental  data  given  by  Schmitt  et  al.  [86].  The 
location  of  z/b  =  0.2  is  near  the  root,  and  z/b  =  0.99  is  at  wing  tip. 
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Overall,  very  good  agreement  is  obtained  between  the  computation  and  experiment 
for  each  cross-section. 


8.8.2  Validation  of  Structural  Solver 

To  validate  the  structural  model  used  in  the  present  study,  the  dynamic  responses 
of  a  flexible  plate  wing  shown  in  Figure  12.135  is  calculated  and  compared  with  the 
results  by  using  the  finite  element  solver  ANSYS.  The  purpose  of  the  study  is  to  find 
out  how  many  mode  shapes  are  required  for  accurate  representation  of  the  structural 
motion  under  dynamic  force. 

The  plate  wing  has  the  same  outline  as  the  AGARD  wing  445.6,  and  its  first 
mode  natural  vibration  frequency  is  nearly  the  same  as  the  corresponding  one  of  the 
AGARD  wing  445.6.  The  thickness  of  the  plate  is  0.3”,  the  root  chord  is  21.96”,  the 
tip  chord  is  14.49”,  and  the  spanwise  length  is  30”.  The  plate  wing  consists  of  80 
elements  and  861  node  points  on  each  side  of  the  wing.  The  plate  wing  is  held  fixed 
on  the  root. 

A  time-dependent  force  is  exerted  at  node  point  510  which  is  located  at  the  mid¬ 
point  of  the  wing  tip.  The  three  components  of  the  force  in  the  unit  of  pound  are: 


fx  =  0.5sm(27r/et),  fy  =  0.3sin(2nfet),  fz  =  0.8sm(27r/et)  (8.22) 

where  the  exciting  frequency  fe  is  equal  to  10  Hz.  The  modal  damping  ratio  Q  =  0.01, 
the  time  step  used  is  0.0005  second.  The  dynamic  responses  at  several  locations  are 
recorded.  Figure  12.136  shows  the  time  histories  of  the  responses  at  the  node  point 
491  which  is  located  almost  at  the  center  of  the  plate  wing.  The  numerical  predictions 
by  the  present  structural  solver  with  first  five  mode  shapes  agree  excellently  with  the 
results  using  ANSYS  with  first  five  mode  shapes  and  the  full  model.  The  three  results 
are  virtually  duplicated  to  each  other  and  are  indistinguishable. 


8.8.3  AGARD  Wing  445.6  Flutter 

The  AGARD  445.6  wing  is  selected  to  demonstrate  the  capability  of  the  present  solver 
for  predicting  the  flutter  boundary.  This  wing  has  a  quarter-chord  sweep  angle  of  45 
degree,  an  aspect  ratio  of  1.65,  a  taper  ratio  of  0.6576,  and  a  NACA65A004  airfoil 
section  in  the  streamwise  direction.  The  weakened  wing  model  (Model  3)  listed  in 
[55]  is  chosen  for  this  study.  The  geometry  of  the  wing  and  its  first  six  mode  shapes 
as  well  as  the  experimental  flutter  results  are  also  provided  in  the  same  report  [55]. 
The  wing  structure  is  modeled  by  its  first  five  natural  vibration  modes  in  the  present 
computation. 
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The  simulations  start  with  the  stationary  rigid  body  wing  model.  After  the  steady 
state  flow  field  around  the  wing  is  fully  developed,  the  rigid  body  wing  is  switched  to 
the  flexible  wing  model.  As  a  small  imposed  perturbation,  the  first  mode  displacement 
of  the  structural  motion  is  set  into  sinusoidal  motion  for  one  cycle  with  the  maximum 
amplitude  of  0.0005  -  0.001  and  the  first  mode  frequency  of  the  wing.  Then  the  wing 
is  allowed  to  deflect  in  response  to  the  dynamic  force  load.  Within  each  physical  time 
step,  the  solution  is  usually  converged  within  50  iterations. 

In  Figures  12.137  through  12.139  the  time  histories  of  generalized  displacements 
of  the  AGARD  wing  445.6  at  =  0.96  are  plotted  for  three  different  V*.  In  these 
figures,  from  V*  =  0.26  to  V*  =  0.315,  the  plots  correspond  to  the  damped,  neutral, 
and  diverging  responses,  respectively.  When  the  value  of  V*  is  smaller  than  the 
critical  value  on  the  flutter  boundary,  the  amplitudes  of  all  modes  decrease  in  time 
corresponding  to  the  damped  response  as  shown  in  Figure  12.137.  Once  the  value 
of  V*  coincides  with  or  is  close  to  the  critical  value,  the  neutral  response  appears  as 
shown  in  Fig.  12.138.  When  the  value  of  V*  is  above  the  neutral  stability  point, 
the  amplitudes  of  first  five  modes  grow  very  fast,  a  diverging  response  is  reached  as 
shown  in  Fig.  12.139. 

For  a  given  Mach  number,  several  runs  with  different  V*  are  needed  to  determine 
the  location  of  the  flutter  boundary  using  bisection  method.  When  V*  is  varied, 
the  free-stream  Reynolds  number  is  changed  accordingly.  Strictly  speaking,  the  free- 
stream  Reynolds  number  needs  to  be  updated  and  the  initial  steady-state  flow  field 
with  actual  Reynolds  number  should  be  re-generated  for  each  run.  In  present  sim¬ 
ulation,  the  initial  flow  field  and  the  Reynolds  number  remain  unchanged  when  V* 
is  varied  since  the  effect  on  final  solution  due  to  small  variation  in  the  free-stream 
Reynolds  number  is  negligible  when  a  turbulence  model  is  used. 

The  comparison  of  computed  flutter  boundary  and  experimental  data  for  AGARD 
Wing  445.6  is  illustrated  in  Figure  12.140.  Overall,  the  computed  results  are  in  good 
agreement  with  the  experimental  data.  The  “sonic  dip”  near  Mach  =1.0  measured  in 
the  experiment  is  very  well  captured  by  the  computation.  The  discrepancy  between 
computed  results  and  experimental  data  may  be  due  to  the  inadequacy  of  the  turbu¬ 
lent  modeling  to  capture  the  shock/wave  boundary  layer  interaction  or  may  be  due 
to  the  inaccurate  measurement  in  the  experiment  as  suspected  by  some  researchers. 
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Chapter  10 
Conclusions 


In  this  research,  a  new  high  efficiency  low  diffusion  upwind  scheme  is  developed  and  is 
applied  to  moving  grid  systems.  The  fully  coupled  fluid-structural  interaction  with  the 
new  upwind  scheme  is  successfully  developed  and  predicts  the  3d  AGARD  wing  flutter 
boundary  very  well.  The  SNM  model  at  time  domain  is  developed  and  validated.  The 
forced  vibration  of  the  NASA  flutter  cascade  is  simulated  with  reasonable  agreement 
with  the  experiment.  A  set  of  nonlinear  non-reflective  boundary  conditions  for  Navier- 
Stokes  equations  on  the  generalized  coordinates  is  developed  to  avoid  the  reflective 
waves.  The  code  is  ready  to  calculate  the  fluid-structural  interaction  of  the  mistunned 
rotor  bladed  disk.  However,  since  the  funding  is  only  3  years,  which  is  one  year  shorter 
than  the  4  years  time  period  originally  proposed,  the  mistunned  rotor  simulation  is 
not  completed  and  will  be  done  in  future  when  the  funding  is  available. 

The  following  is  the  specific  conclusions  achieved  in  this  research: 


10.1  The  New  E-CUSP  scheme 

A  new  efficient  upwind  scheme  based  on  the  concept  of  convective  upwind  and  split 
pressure  (CUSP)  is  developed.  The  upwinding  of  the  convective  term  and  the  pressure 
splitting  are  consistent  with  their  characteristics  directions.  The  numerical  dissipation 
of  the  new  scheme  at  stagnation  is  low  and  is  not  greater  than  that  of  the  Roe  scheme. 
The  scheme  hence  is  able  to  resolve  accurately  wall  boundary  layers,  and  are  able  to 
capture  crisp  shock  waves  and  exact  contact  discontinuities.  The  performance  of  the 
new  scheme  is  compared  with  the  Roe  scheme,  AUSM+  scheme,  Van  Leer  scheme, 
and  Van  Leer-Hanel  scheme. 

For  the  ID  Sod  shock  tube  problem  using  Euler  explicit  scheme,  the  new  scheme 
has  the  crispest  shock  profile  and  highest  allowable  CFL  number  of  1.0.  For  a  slowly 
moving  contact  surface,  the  new  scheme  is  demonstrated  to  capture  the  exact  contact 
surface  discontinuity  with  the  maximum  allowable  CFL  of  1.0,  which  is  far  greater 


87 


88 


CHAPTER  10.  CONCLUSIONS 


than  that  of  the  other  schemes.  For  a  quasi- ID  transonic  nozzle,  all  the  other  schemes 
generate  expansion  shocks  at  the  sonic  point.  The  new  scheme  does  not  have  the 
expansion  shock  even  though  it  has  a  glitch  at  the  sonic  point,  which  is  due  to  the 
discontinuity  of  the  first  derivative  of  the  pressure  splitting  at  sonic  point. 

For  a  Mach=2.0  supersonic  adiabatic  laminar  flat  plate  boundary  layer,  the  new 
scheme  is  able  to  accurately  resolve  the  boundary  layer  velocity  and  temperature 
profiles  using  the  first  order  differencing.  The  solution  is  as  accurate  as  that  of  the 
Roe  scheme  and  the  AUSM+  scheme  and  hence  demonstrates  the  low  diffusion  of  the 
new  scheme. 

For  a  transonic  converging-diverging  nozzle,  oblique  shock  waves  and  reflections 
are  crisply  captured  even  though  the  shock  waves  do  not  align  with  the  mesh  lines. 
The  predicted  wall  surface  isentropic  Mach  number  distribution  agrees  well  with  the 
experiment.  For  a  transonic  inlet-diffuser  with  shock/turbulent  boundary  layer  inter¬ 
action,  the  new  scheme  and  the  Roe  scheme  predict  the  surface  pressure  distributions 
agreeing  well  with  the  experiment  for  the  case  of  a  weak  shock.  For  the  strong  shock 
case,  both  the  new  scheme  and  the  Roe  scheme  over  predict  the  strength  of  the  shock 
wave.  However,  the  pressure  distribution  predicted  by  the  new  scheme  is  closer  to 
the  experiment.  The  AUSM+  solution  has  large  pressure  oscillations. 

In  the  applications  of  the  original  Zha-Hu  E-CUSP  scheme,  it  is  found  that  there  is 
temperature  oscillations  near  wall  boundaries.  The  scheme  is  modified  by  replaceing 
the  pressure  term  with  the  total  enthalpy  in  the  smoothing  factor  of  the  energy 
equation.  The  temperature  oscillations  are  removed  in  the  modified  scheme,  which  is 
named  Zha  E-CUSP2  scheme. 

10.2  Non- Reflective  Boundary  Conditions 

The  non-reflective  boundary  conditions  of  Poinsot  and  Lele[33]  for  3D  Navier-Stokes 
equations  are  extended  to  the  generalized  coordinates  in  this  paper.  The  character¬ 
istic  form  Navier-Stokes  equations  in  conservative  variables  are  given.  The  NRBC  is 
applied  numerically  in  an  implicit  time  marching  method.  The  governing  equations 
for  inner  domain  and  NRBC  are  solved  simultaneously  in  a  fully  couple  manner. 

For  the  unsteady  subsonic  vortex  propagating  flow,  the  fixed  pressure  outflow 
boundary  conditions  imposing  the  exit  pressure  generates  serious  wave  reflection  and 
the  flow  field  is  distorted,  whereas,  the  NRBC  developed  in  this  paper  generates  clean 
results  with  no  wave  reflection  and  solution  distortion. 

For  the  transonic  inlet-diffuser,  the  NRBC  is  not  necessary  for  steady  state  solu¬ 
tions  since  the  reflective  waves  are  diffused  when  the  solutions  are  converged.  How¬ 
ever,  for  unsteady  flows,  the  NRBC  is  essential.  The  fixed  pressure  outflow  boundary 
conditions  generate  strong  reflective  waves  due  to  the  shock/boundary  layer  inter- 
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action,  which  makes  the  shock  oscillating  motion  far  greater  than  the  experimental 
data.  When  the  NRBC  is  applied,  the  shock  oscillation  is  dramatically  reduced  and 
the  computed  time  averaged  pressure  distributions  and  frequency  spectrum  agree 
much  better  with  the  experiment  than  FPBC. 

The  NRBC  is  essential  for  unsteady  fluid-structural  interactions,  in  particular  for 
the  internal  flows  such  as  the  mistunned  rotors. 


10.3  2D  Flutter  Prediction 

The  efficient  high  resolution  E-CUSP  upwind  schemes  developed  in  this  research 
have  been  successfully  extended  and  applied  to  calculate  flow-induced  vibration  with 
a  moving  grid  based  on  fully  coupled  fluid-structural  interaction. 

For  an  elastically  mounted  cylinder,  various  cases  with  different  structural  param¬ 
eters  have  been  calculated.  The  predicted  displacement  agrees  very  well  with  the 
experiment  and  the  numerical  results  of  other  researchers. 

For  the  forced  pitching  NACA  64A010  airfoil,  the  computed  lift  oscillation  agrees 
accurately  with  the  experiment.  The  computed  moment  oscillation  has  large  devia¬ 
tion  from  the  experiment,  but  the  results  have  a  similar  order  of  accuracy  as  other 
researchers  have  achieved.  The  discrepancy  may  be  due  to  the  inaccurate  prediction 
of  the  surface  shear  stress  caused  by  the  inadequacy  of  the  turbulence  modeling  and 
experimental  uncertainties. 

The  same  airfoil  has  been  calculated  as  an  elastically  mounted  airfoil  with  free 
vibration  in  plunging  and  pitching  directions.  The  computations  have  been  carried 
out  using  different  values  of  velocity  index  V*  which  are  below,  close,  and  beyond 
the  critical  point  on  the  flutter  boundary  when  Mach  number  is  fixed.  The  corre¬ 
sponding  responses  of  the  airfoil  flows  are  well  simulated.  The  predicted  value  of  V* 
at  the  bottom  of  the  transonic  dip  is  consistent  with  the  numerical  results  of  other 
researchers  [5,  73,  70], 


10.4  SNM  model  of  Transient  Response 

The  transient  response  (time  domain)  structural  vibration  model  for  mistunned  rotor 
bladed  disk  based  on  the  efficient  Subsets  of  Nominal  Modes  (SNM)  model  is  devel¬ 
oped.  The  vibration  response  results  predicted  by  the  SNM  model  for  a  full  annulus 
bladed  disk  with  blade  frequency  variation  agree  well  with  the  results  predicted  by  the 
finite  element  model.  The  simplified  model  was  designed  such  that  transient  response 
computations  can  be  performed  in  ANSYS  with  tractable  computational  cost.  The 
comparison  of  the  results  showed  that  SNM  model  can  capture  both  low  frequency 
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and  high  frequency  responses.  In  the  case  of  low  frequency  excitation,  only  one  family 
of  blade  modes  is  needed  in  SNM.  In  the  case  of  high  frequency  excitation,  multiple 
families  of  blade  modes  are  needed  in  SNM.  The  computational  efficiency  was  im¬ 
proved  by  at  least  30  times  over  ANSYS  direct  simulations  for  the  cases  examined. 
The  improvement  is  expected  to  be  much  greater  when  high  fidelity  ANSYS  models 
are  needed. 


10.5  Separated  Flows  in  NASA  3D  Flutter  Cas¬ 
cade 

1.  The  high  incidence  cascade  separation  flow  shows  a  sinusoidal  pattern  on  the 
oscillation  of  the  surface  pressure  and  the  separation  bubble  size. 

2.  The  leading  edge  vortex  shedding  is  the  mechanism  behind  the  unsteady  char¬ 
acteristics  of  the  subsonic  high  incidence  separation  flow.  New  vortexes  are 
continuously  generated  at  the  suction  surface  leading  edge.  The  new  vortex 
grows  and  pushes  the  old  vortexes  downstream.  The  interaction  between  the 
vortexes  results  in  the  periodical  oscillation  of  the  separation  bubble  size  and 
the  surface  pressure.  The  vortex  generation,  pressure  variation  and  separation 
length  oscillation  have  the  same  frequency  characteristics  with  a  phase  differ¬ 
ence. 

3.  The  characteristics  of  the  separation  flow  is  determined  by  the  inlet  Mach  num¬ 
ber.  When  the  inlet  flow  goes  from  lower  subsonic  to  higher  subsonic,  the 
size  and  the  oscillation  intensity  of  the  separation  bubble  are  enhanced.  The 
flowfield  oscillation  peak  frequency  increases.  When  the  inflow  goes  further  to 
supersonic,  the  flow  is  attached  on  the  leading  edge.  A  small  size  separation 
bubble  due  to  the  interaction  of  the  shock  wave  and  the  turbulent  boundary 
layer  is  located  right  after  the  shock  wave. 


10.6  Forced  Vibration  of  the  NASA  Flutter  Cas¬ 
cade 

The  fully  nonlinear  time  dependent  Navier-Stokes  equations  are  solved  with  the  paral¬ 
lel  computation  technique  to  simulate  the  unsteady  flow  field  in  a  full  scale  compressor 
cascade  with  forced  blade  vibration.  The  calculation  in  this  paper  is  conducted  with 
a  low  incidence  of  0°  and  a  subsonic  inflow  M—  0.5.  The  blade  motion  amplitude 
is  1.2°and  the  inter  blade  phase  angle  is  180°.  The  full  scale  computation  is  carried 
out  for  3  reduced  frequencies,  0.4,  0.8  and  1.2.  The  blade  stability  under  different 
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vibration  frequency  is  analyzed.  The  end  wall  influence  on  the  steady  and  unsteady 
flow  characteristics  is  studied  by  comparing  the  full  scale  results  with  the  2-passage 
periodic  cascade  results  at  kc  =0.8.  The  conclusions  are  the  following: 

1.  The  flow  pattern  in  the  full  scale  cascade  shows  that  the  fiowfield  is  affected  by 
the  existence  of  the  end  walls.  The  steady  state  blade  surface  pressure  varies 
with  the  the  blade  position  in  the  cascade.  The  periodicity  of  the  flow  pattern 
is  improved  by  adjusting  the  end  wall  configuration.  The  end  wall  influence 
attenuates  rapidly  from  boundary  passages  to  center  passages.  Good  periodicity 
is  achieved  in  the  inner  passages.  The  full  scale  computation  gives  better  results 
of  the  unsteady  pressure  coefficient,  local  stability  and  aerodynamic  moment 
among  the  center  blades. 

2.  All  blades  in  the  full  scale  cascade  are  stable,  which  is  indicated  by  a  positive 
damping  coefficient.  The  damping  coefficient  is  more  uniformly  distributed 
on  center  blades.  The  most  stable  and  the  least  stable  blades  are  the  two 
boundary  blades.  The  damping  coefficient  and  its  variation  across  the  center 
blades  increase  with  the  increasing  vibration  frequency. 

3.  The  unsteady  pressure  coefficients  are  predicted  well  compared  with  the  ex¬ 
periment  measurement.  The  local  stability  trend  is  correctly  predicted  in  the 
numerical  computation.  The  blade  local  stability  is  over-predicted  on  LE  and 
under-predicted  on  the  aft  part.  The  destabilization  region  located  at  the  front 
part  of  the  blade  is  predicted  smaller  compared  with  the  experiment.  The  pre¬ 
dicted  stabilization  and  destabilization  increase  with  the  increasing  frequency, 
so  do  the  damping  coefficients. 


10.7  3D  AGARD  Wing  Flutter  Prediction 

A  numerical  methodology  with  fully  coupled  fluid-structural  interaction  for  predi¬ 
cating  3-D  transonic  wing  flutter  has  been  developed.  A  dual-time  step  implicit 
unfactored  Gauss-Seidel  iteration  with  Roe  scheme  are  employed  in  the  flow  solver. 
A  modal  approach  structure  solver  is  used  to  simulate  the  wing  response.  An  effi¬ 
cient  mesh  deformation  strategy  based  on  an  algebraic  method  is  developed  and  is 
shown  to  be  accurate  and  robust.  The  flow  and  structure  solvers  are  fully  coupled  vis 
successive  iterations  within  each  physical  time  step. 

The  accuracy  of  the  modal  approach  structure  solver  has  been  verified  by  using 
the  finite  element  solver  ANSYS.  The  results  indicate  that  the  first  five  modes  are 
sufficient  to  accurately  model  the  wing  structure  in  the  present  study. 

The  computed  flutter  boundary  of  AGARD  wing  445.6  for  free  stream  Mach  num¬ 
bers  ranging  from  0.499  to  1.141  is  presented  and  compared  well  with  the  experimental 
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data  except  for  =  1.171,  where  the  flutter  boundary  is  over-predicted.  The  sonic 
dip  is  very  well  captured. 


10.8  Future  Work 

With  the  foundation  laid  in  the  research,  the  next  step  is  to  apply  the  fully  coupled 
fluid-structural  interaction  solver  to  mistunned  rotor  bladed  disk. 
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Figure  12.1:  Temperature  distribution  of 
the  Sod  ID  shock  tube  computed  by  Zha 
E-CUSP  scheme 


Figure  12.3:  Temperature  distribution  of 
the  Sod  ID  shock  tube  computed  by  Van 
Leer  scheme 


Figure  12.2:  Temperature  distribution  of  Figure  12.4:  Temperature  distribution  of 
the  Sod  ID  shock  tube  computed  by  Roe  the  Sod  ID  shock  tube  computed  by  Van 
scheme  Leer-Hanel  scheme 
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Figure  12.5:  Temperature  distribution 
of  the  Sod  ID  shock  tube  computed  by 


AUSM+  scheme 


Figure  12.7:  Computed  density  profile  of 
a  slowly  moving  contact  surface  using  the 
Roe  scheme 


Figure  12.6:  Computed  density  and  ve-  Figure  12.8:  Computed  density  and  ve¬ 
locity  profiles  of  a  slowly  moving  contact  locity  profiles  of  a  slowly  moving  contact 
surface  surface 
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Figure  12.10:  Computed  velocity  profiles 
of  the  laminar  boundary  layer  using  1st 
order  schemes 


Figure  12.12:  Computed  Mach  number 
contours  using  the  Zha  CUSP  scheme 
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Figure  12.15:  Static  pressure  distribution 
Figure  12.13.  Adiabatic  Mach  number  computed  on  the  upper  surface  of  the 
distribution  computed  on  the  wall  surface  inlet-diffuser  (pout/pt  —  -83) 
of  the  nozzle 


Figure  12.14:  Computed  Mach  number 
contours  using  the  Zha  CUSP  scheme 
with  Pout / Pt  =  .83 


Figure  12.16:  Static  pressure  distribution 
computed  on  the  upper  surface  of  the 
inlet-diffuser  ( Pout/Pt  =  -72) 
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Figure  12.17:  Comparison  of  computed 
pressure  contours  using  the  Zha  CUSP 
scheme,  Roe  Scheme,  and  Liou’s  AUSM+ 
scheme  ( Pout/Pt  =  -72) 
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Figure  12.19:  Zoomed  temperature  con¬ 
tours  near  the  symmetric  boundary  us¬ 
ing  the  Zha  CUSP(a)  and  Zha  CUSP2 
scheme(b) 


Figure  12.18:  Computed  Mach  number 
contours  using  the  Zha  CUSP  (a)  and  Zha 
CUSP2  scheme  (b) 
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Figure  12.20:  Sketch  of  the  elastically  F[^re  12  22:  The  near  wal1  zone  mesh 
mounted  cylinder  around  the  solid  surface  of  the  cylinder. 


Figure  12.21:  Sketch  of  the  elastically 
mounted  airfoil. 


Figure  12.23:  Time  histories  of  the  lift 
and  drag  coefficients  of  the  stationary 
cylinder  due  to  vortex  shedding. 
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Iteration  number 

Figure  12.24:  Typical  convergence  his¬ 
tories  of  the  CFD  and  structural  solvers 
within  one  physical  time  step. 
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Figure  12.26:  The  trajectory  of  the  center 
position  of  the  oscillating  cylinder,  /is  = 
1.2732,  C  =  0.63326. 


-5 Log10  of  reduced  damping 
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Figure  12.27:  Comparison  of  the  com- 
Figure  12.25:  Vorticity  contours  around  puted  amplitude  with  Griffin’s  experi- 
oscillating  cylinder,  fis  =  1.2732,  £  =  mental  data  for  the  elastically  mounted 
0.03166.  cylinder. 
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Figure  12.28:  Surface  pressure  coeffi¬ 
cient  comparison  for  RAE  2822  airfoil  at 
M=0.729. 
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Figure  12.30:  Comparison  of  computed 
lift  coefficient  with  Davis’  experimental 
data  for  the  forced  pitching  airfoil. 


Figure  12.31:  Comparison  of  computed 
Figure  12.29:  The  near  wall  O-type  mesh  moment  coefficient  with  Davis’  experi- 
around  the  NACA  64A010  airfoil.  mental  data  for  the  forced  pitching  airfoil. 
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Figure  12.32:  Comparison  of  the  mean  Figure  12-34:  Time  histories  of  plung- 
pressure  distribution  along  the  pitching  *ng  and  pitching  displacements  for  = 
ajrfoij  0.825  and  V*  —  0.55  -  Damped  response. 


Figure  12.33:  Comparison  of  the  instan¬ 
taneous  pressure  distributions  at  various 
angles  of  attack. 


Figure  12.35:  Time  histories  of  plung¬ 
ing  and  pitching  displacements  for  Mm  = 
0.825  and  V*  =  0.615  -  Neutrally  stable 
response. 
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Figure  12.36:  Time  histories  of  plung¬ 
ing  and  pitching  displacements  for  Mm  = 
0.825  and  V*  =  0.70  -  Diverging  response. 


Figure  12.38:  Comparison  of  computed 
flutter  boundaries  -  — -  versus  Mach  num- 

“o 

ber. 
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Figure  12.37:  Comparison  of  computed 
flutter  boundaries  -  Speed  index  versus 
Mach  number. 


Figure  12.39:  Time  histories  of  plung¬ 
ing  and  pitching  displacements  for  = 
0.925  and  V*  =  5.5  -  Limit  Cycle  Oscilla¬ 
tion. 
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...  .  ,  ,  .  Figure  12.42:  Finite  element  model  of  24- 

Figure  1240:  Time  histories  of  plunging  ,  °  ,. 

and  pitching  displacements  for  =  0.9  a  6  1S 
and  V*  =  2.5  with  the  second  mode  vi¬ 
bration  dominant. 


Dimensionless  Structural  Time 


(Prescribed)  Force  at  Blade  Tip  (0  Nodal  Diameter) 
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Figure  12.41:  Time  histories  of  plung¬ 
ing  and  pitching  displacements  for  = 
0.875  and  V*  =  2.5  -  ‘Standing’  status. 


Figure  12.43:  Prescribed  force  at  blade 
tip,  zeroO  nodal  diameter  at  low  fre¬ 
quency. 
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r>  ,  ,  .  , ,  ,  Figure  12.47:  Comparison  of  the  pre¬ 
figure  12.45:  Prescribed  force  at  blade  ,  ,  ,  .  .  ,  .  . 

, .  ,  ,  ,  ,  ,  ,  dieted  transient  response  at  low  ire- 

tip,  zero  nodal  diameter  at  low  frequency. 

quency. 
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Response  of  Blade  1  Tip  (Prescribed)  Force  a(  Blade  Tip  (0  Nodal  Diameter) 


Figure  12.48:  Comparison  of  the  pre-  Figure  12.50:  Prescribed  force  at  blade 
dieted  transient  response  at  low  fre-  tip,  zero  nodal  diameter  at  high  fre¬ 
quency.  quency. 


Figure  12.49:  Prescribed  force  at  blade 
tip,  zero  nodal  diameter  at  high  fre¬ 
quency. 


Figure  12.51:  Prescribed  force  at  blade 
tip,  zero  nodal  diameter  at  high  fre¬ 
quency. 
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Figure  12.52:  Comparison  of  the  pre-  Figure  12.54:  Comparison  of  the  pre¬ 
dicted  transient  response  at  high  fre-  dieted  transient  response  at  high  fre¬ 
quency.  quency. 
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Figure  12.53:  Comparison  of  the  pre-  Figure  12.55:  Comparison  of  the  pre¬ 
dicted  transient  response  at  high  fre-  dieted  transient  response  at  high  fre¬ 
quency.  quency. 
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t  =  0.5  t  =  1 .0 


t  =  1 .5  t  =  2.0 


Figure  12.56:  Comparison  of  the  pre¬ 
dicted  transient  response  at  high  fre¬ 
quency. 


Figure  12.58:  Density  contours  at  four  in¬ 
stants  for  a  vortex  leaving  the  domain  us¬ 
ing  imposed  exit  static  pressure  boundary 
conditions. 


t*  =  0.5  t"  =  1.0 


t*  *  1.5  t’  =  2.0 


Figure  12.57:  Comparison  of  the  pre¬ 
dicted  transient  response  at  high  fre-  Figure  12.59:  Density  contours  at  four  in- 
quency.  stants  for  a  vortex  leaving  the  domain  us¬ 

ing  NRBC  exit  boundary  conditions. 
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t  =  0.5  t‘=  1.0 


Figure  12.61:  ( uT  —  uT00)/uT00  contours 
at  three  instants  for  a  vortex  leaving  the 
domain  using  NRBC  exit  boundary  con¬ 
ditions. 
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Figure  12.67:  Steady  state  static  pres¬ 
sure  distribution  along  the  bottom  wall 
for  Poutut/Pt  =  0.72. 


Figure  12.69:  Time  averaged  unsteady 
pressure  distribution  along  the  bottom 
wall  for  Pvutiet/Pt  ~  0.72. 


Figure  12.68:  Time  averaged  unsteady  Figure  12.70:  Power  spectral  density  for 
pressure  distribution  along  the  top  wall  the  static  pressure  fluctuations  at  x/H  = 
for  Poutiet/Pt  =  0.72.  14.218,  p0utiet/Pt  =  0.72. 
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(a)  bottom  (b)  midspan  (c)  top 
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Figure  12.73:  Bottom,  mid-span  and  top 
planes  of  the  cascade  mesh 
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Figure  12.72:  Cascade  3D  mesh 


Figure  12.75:  Mid-span  Stream  Lines  at 
Incidence  Angle  10° 
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Figure  12.76:  Mid-span  Suction  Surface 
Streamlines  at  Incidence  Angle  10° 
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Figure  12.79:  Separation  Zone  Length 
Variation  Frequency  Spectrum 


Figure  12.77:  Mid-span  Surface  Pressure 
Distribution  at  Incidence  Angle  10° 


Figure  12.80:  Suction  Surface  Check 
Point  Pressure  Variation  with  Time 


Figure  12.78:  Separation  Zone  Length 
Variation  with  Time 
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a,  t  =4.4305  ms,  8.4tc 
8.8  tc 


Figure  12.81:  Suction  Surface  Check 
Point  Pressure  Frequency  Spectrum 


t  =4.6415  ms, 


0.76 

0.755 

0.75 

0.745 

0.74 

0.735 

0.73 

0.725 


Figure  12.82:  Check  Point  Pressure  and 
Separation  Locations  Relation 


e,  t  =5.2744  ms,  10.0tc  f,  t  =5.4854  ms, 
10.4fc 


g,  t  =5.6964  ms,  10.8fc  h,  t  =5.9073  ms, 
11.2ic 


Figure  12.83:  Separation  Bubble  Evolu¬ 
tion 
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Figure  12.84:  Separation  Zone  Length  Os-  Figure  12.87:  Checkpoint  Pressure  Oscil- 
cillation  (Ma= 0.8)  lation  Frequency  Spectrum  (Ma=0.8) 


Figure  12.85:  Separation  Zone  Length  Os-  Figure  12.88:  Checkpoint  Pressure  Oscil- 
cillation  Frequency  Spectrum  (Ma=0.8)  lation  (Ma— 1.18) 


Figure  12.86:  Checkpoint  Pressure  Oscil-  Figure  12.89:  Checkpoint  Pressure  Oscil¬ 
lation  (Ma= 0.8)  lation  Frequency  Spectrum  (Mo=1.18) 
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Figure  12.90:  Cascade  structure  and  CFD 
setup 


Blade  Incidence  (degree) 


Pressure  Oscillation  on  Lower  Surface 


Figure  12.91:  Cascade  Mesh 


Pressure  Oscillation  on  Upper  Surface 


Figure  12.94:  Pressure  oscillation  at  dif¬ 
ferent  chordwise  locations 


Figure  12.92:  Additional  mesh  layer  for 
better  orthogonality 


Figure  12.93:  Pressure  oscillation  history 


Figure  12.95:  Mach  number  contours  of 
the  cascade 
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Figure  12.96:  Mach  number  contours  at 
different  time  steps 


Figure  12.97:  Mach  number  contours  at 
different  time  steps 


Figure  12.98:  Mach  number  contours  at 
different  time  steps 


Figure  12.99:  Mach  number  contours  at 
different  time  steps 


Figure  12.100:  Mach  number  contours  at 
different  time  steps 


Figure  12.101:  Mach  number  contours  at 
different  time  steps 
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Figure  12.103:  Mach  number  contours  at 
different  time  steps 


Figure  12.104:  Mach  number  contours  at 
different  time  steps 


Figure  12.105:  Mach  number  contours  at 
different  time  steps 
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Figure  12.106:  Unsteady  pressure  analy¬ 
sis 
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Figure  12.107:  Computation  Domain 
Configuration 


Figure  12.108:  Computation  mesh 


Figure  12.109:  Two-passage  cascade  peri¬ 
odic  computation  domain 
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Figure  12.110:  Steady  state  pressure  co¬ 
efficient  in  2-passage  cascade 
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Figure  12.111:  Pressure  oscillation  his¬ 
tory,  kc=  0.8 


Figure  12.112:  Unsteady  pressure  coeffi¬ 
cient  on  upper  surface  in  2-passage  cas¬ 
cade,  kc  =  0.8 
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Figure  12.113:  Unsteady  pressure  coeffi-  FiSure  ®  sta‘e  Pressurc  c°- 

cient  on  lower  surface  in  2-passage  cas-  e  c'eri*  m  u  sca  e  casca  e 
cade,  kc=  0.8 


Figure  12.117:  Unsteady  state  pressure 
coefficient  on  upper  surface,  kc  =0.8 


Figure  12.114:  Local  stability  analysis  in 
2-passage  cascade,  fcc=0.8 


x/C 


Figure  12.115:  Steady  state  Mach  con-  Figure  12.118:  Unsteady  state  pressure 
tours  in  full  scale  cascade  coefficient  on  lower  surface,  kc  =0.8 
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Figure  12.124:  Unsteady  pressure  coeffi-  Figure  12.127:  Local  stability  analysis, 
cient  on  lower  cascade,  kc  =1.2  kc  =0.4 


Figure  12.125:  Local  stability  analysis,  Figure  12.128:  Damping  coefficient  distri- 
kc  =1.2  bution,  kc  =0.4 


Figure  12.126:  Damping  coefficient  distri-  Figure  12.129:  Unsteady  aerodynamic 
bution,  kc  =1.2  moment  oscillation  comparison 
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Figure  12.130:  Local  stability  analysis 
comparison 


Figure  12.132:  Fully  coupled  flow- 
structure  interaction  calculation  proce¬ 
dure 


Figure  12.131:  A  sketch  of  the  mesh  d 
formation 


Figure  12.133:  The  mesh  around  the  ON- 
ERA  M6  wing  surface. 
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Figure  12.136:  Histories  of  the  dynamic 
Figure  12.134:  Pressure  coefficients  on  the  reSponses  at  node  point  491. 
wing  surface  at  different  cross-section. 


Figure  12.135:  Plate  wing  geometry. 


Dimensionless  time 


Figure  12.137:  Time  histories  of  the  gen¬ 
eralized  displacements  of  first  five  modes 
for  Moo  =  0.96  and  V*  =  0.26  -  Damped 
response. 
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Figure  12.138:  Time  histories  of  the  ge 
eralized  displacements  of  first  five  mod 
for  Mqo  =  0.96  and  V*  =  0.285  -  Neutral 
stable  response. 


Figure  12.139:  Time  histories  of  the  gen¬ 
eralized  displacements  of  first  five  modes 
for  Mqo  =  0.96  and  V*  =  0.315  -  Diverging 


response. 
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Figure  12.140:  Comparison  of  predicted 
wing  flutter  boundary  with  experimental 
data  for  AGARD  Wing  445.6. 


